- Email: [email protected]
Sub-grid scale model for convection-driven dynamos in a rotating plane layer
Physics of the Earth and Planetary Interiors 153 (2005) 108–123
Sub-grid scale model for convection-driven dynamos in a rotating plane layer Hiroaki Matsui∗ , Bruce A. Buffett Department of Geophysical Sciences, The University of Chicago, 5734 S Ellis Avenue, Chicago, IL 60637, USA Received 2 November 2004; accepted 3 March 2005
Abstract We develop a sub-grid scale (SGS) model for convection-driven dynamos, based on the nonlinear gradient model of Leonard [Leonard, A., 1974. Energy cascade in large-eddy simulations of turbulent fluid flows. Adv. Geophys. 18, 237–248]. The predictions of the SGS model are tested using snapshots from a direct numerical simulation (DNS). Good agreement is obtained throughout most of the fluid, but discrepancies occur near the boundaries. We also implement the SGS model in a large-eddy simulation (LES) to assess the temporal behavior of the SGS model. Both the DNS and LES have large time variation in the kinetic and magnetic energies, so comparisons are made using time averages. The time-averaged energies in the LES are slightly larger than those in the resolved DNS, but not as high as those obtained in an unresolved DNS with the same grid resolutions as the LES. The source of dissipation in the LES is revealed by examining the energy balances. The work done by the SGS Lorentz term and the magnetic energy generated by the SGS induction term are both negative, indicating that these SGS terms transfer energy from the resolved scales into the unresolved scales. © 2005 Elsevier B.V. All rights reserved. Keywords: Geodynamo; Magnetohydrodynamics; Sub-grid scale model
1. Introduction In the last 10 years, geodynamo simulations have begun to reveal the basic processes responsible for generating the geomagnetic field. Following the initial studies of Glatzmaier and Roberts (1995) and Kageyama and Sato (1995), a number of additional geodynamo models have been developed, which reproduce many characteristics of the geomagnetic field (e.g. Kuang and Bloxham, 1997, 1999; Christensen et al., 1999; Sakuraba and Kono, 1999). However, the fluid motions in the outer core are expected to have a broad range of length scales, ranging from the dissipative scale in the boundary layer (less than L ∼ 10 m) to the geometry of the outer core ∗
Corresponding author. Tel.: +1 773702 8107; fax: +1 773702 9505. E-mail address: [email protected] (H. Matsui)
0031-9201/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2005.03.019
(L ∼ 106 m). Current numerical simulations cannot use realistic parameters for the Earth’s core because of limited spatial resolution. To account for the influence of unresolved flow, all simulations have so far used turbulent diffusivities (e.g. Glatzmaier and Roberts, 1995; Sakuraba and Kono, 1999). Some of these studies use constant diffusivities, whereas others allow the diffusivities to depend on the length scale of the motion (e.g. hyperdiffusivity). Zhang and Jones (1997) noted that hyperdiffusivity alters the structure of large scale convection, which means that the choice of sub-grid scale models can affect the resolved fields. Small scale convection in the outer core is expected to be strongly anisotropic due to the influences of the Earth’s rotation and magnetic field (Braginsky and Meytlis, 1990; Shimizu and Loper, 1997; Matsushima et al., 1999). However, few models deal with anisotropic
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
diffusivities in numerical calculations. One exception is described by Phillips and Ivers (2001, 2003), who have developed a tensor diffusivity model for use in pseudospectral geodynamo codes. Another example is given by Matsushima (2001), who used magnetoconvection calculation in a small periodic domain of the outer core to devise an anisotropic model for turbulent thermal diffusion. Buffett (2003) also used magnetoconvection calculations to test predictions of SGS heat and momentum fluxes using turbulent viscosity, hyper diffusivity, as well as the Smagorinsky and the similarity models. He concluded that the similarity model successfully reproduces the anisotropy of the SGS heat flux with respect to the rotation axis and direction of the imposed magnetic field. In the present study, we implement a variant of the similarity model in a MHD dynamo model using the finite-element method (FEM). The model is based on the nonlinear gradient model of Leonard (1974), which represents the SGS terms using the gradients of the resolved fields. This approach does not require explicit filtering of the resolved field, as is often the case in SGS modeling (see Lesieur and Metais, 1996; Meneveau and Katz, 2000 for reviews). We apply the nonlinear gradient model to a plane-layer dynamo. Large-eddy simulations (LES) are performed using the nonlinear gradient model and the results are compared with estimates from a direct numerical simulation (DNS) using a finer grid. The goals of the present study are three-fold. First, we attempt to identify and develop practical methods that are capable of reproducing the strong anisotropy expected in the SGS heat and momentum flux. Encouraging results have previously been obtained using the similarity model in magnetoconvection calculations (Buffett, 2003; Matsushima, 2001). However, methods have not yet been tested in dynamo calculations. Second, we need to extend the similarity model to include SGS Lorentz and induction terms. Previous dynamo models in terrestrial (e.g. Glatzmaier and Roberts, 1997) and stellar contexts (e.g. Brun et al., 2004) account for only the eddy diffusivities. The extension of the similarity model to include the SGS Lorentz and induction terms is straightforward, and our calculations show that these additional terms contribute significantly to the evolution of the kinetic and magnetic energy at the resolved scales. In fact, the similarity model can potentially transfer energy into the resolved scales from the unresolved scales. Such a process might emulate the α-effect, which was originally proposed to describe the generation the large scale magnetic field by small scale helical flow (e.g. Parker, 1955; Moffatt and Proctor, 1982; R¨adler et al., 1990). However, the transfer of the energy into the resolved fields may also causes numerical instabilities. Previous stud-
109
ies of homogeneous turbulence have shown that the similarity model does not dissipate enough kinetic energy (e.g. Winckelmans et al., 2001). We test the similarity model in the dynamo problem without the addition of an eddy diffusivity term to address this specific question. We show that the SGS model may actually overestimate the dissipation of magnetic energy in the resolved scales. This and other shortcomings motivate proposals for further improvements to the SGS similarity model. 2. Simulation models 2.1. Sub-grid scale modeling Sub-grid scale modeling often relies on spatial filtering to eliminate field components with length scales smaller than the grid spacing ∆. Each field in the model (e.g. velocity u, . . .) is convolved with a filter function G(x) to define the large scale field (e.g. Leonard, 1974): u¯ = G(x − x )u(x ) dV, (1) where u¯ is called the filtered or resolved velocity. Many types of filter functions are available, and we adopt a Gaussian filter, given by 6 1 6x2 G(x) = exp − 2 . (2) π∆ ∆ The filter is applied to the governing equations to obtain a set of equations for resolved fields. The filtering procedure introduces one new term into the governing equations for each nonlinear term in the original equations (see, for example, Buffett, 2003). We consider the problem of convection-driven motion of the conductive fluid in a plane layer that rotates with the angular velocity = Ωˆz. We use the Boussinesq approximation and take gravity to be a constant vector g = −gˆz. The fluid has a kinetic viscosity ν, thermal diffusivity κ, thermal expansion coefficient α, and magnetic diffusivity η. To obtain the normalized equations, we use the thickness D of the fluid layer as a length scale, and the viscous diffusion time, D2 /ν, as the characteristic time scale. The magnetic field is normalized by B02 = ρµ0 ηΩ, where ρ is the density of the fluid and µ0 is the magnetic permeability. We apply the filter to the basic equations for a conducting fluid with the Boussinesq approximation to obtain ∂u¯ ¯ + ∇(u¯ u) ∂t ¯ + RaE−1 (T¯ − T0 )ˆz = −∇ χ¯ + ∇ 2 u¯ − 2E−1 (ˆz × u) −1 −1 ¯ B)−∇T ¯ + Pm−1 E−1 ∇(B ij +Pm E ∇Mij , (3)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
110
∂T¯ = −∇ · (u¯ T¯ ) + Pr −1 ∇ 2 T¯ − ∇ · I, ∂t
(4)
¯ ∂A ¯ + u¯ × B ¯ + a, = −∇ ϕ¯ + Pm−1 ∇ 2 A ∂t
(5)
¯ = 0, ∇ · u¯ = ∇ · A
(6)
and ¯ ¯ = ∇ × A, B
(7)
where u, χ, T, T0 , B, A, and ϕ are the velocity, modified pressure, temperature perturbation, reference temperature, magnetic field, vector potential of the magnetic field, and arbitrary scalar potential to satisfy the Coulomb gauge, respectively. The bar indicates the filtered fields. The reference temperature T0 is given by T0 = (ztop − z)/D, where z is the vertical coordinate. The Prandtl number Pr, the modified Rayleigh number Ra, the Ekman number E, and the magnetic Prandtl number Pm are defined as follows: ν Pr = , (8) κ αg∆TD , Ων ν E= , ΩD2 Ra =
(9) (10)
and Pm =
ν . η
(11)
Four SGS terms appear in the basic equations (3)–(5), each corresponding to a nonlinear term in the unfiltered equations. These include the SGS heat flux: Ii = ui T − u¯ i T¯ ,
The nonlinear gradient model was original proposed by Leonard (1974). It can now be viewed as a version of the similarity model, which is based on the scale invariance of the turbulence. We describe this model using the SGS heat flux I = uT − u¯ T¯ as an example. The similarity model relies on a second filter with a width that ˜ is greater than the grid size ∆. The second filter G(x) is ¯ applied to the resolved field u(x) to give ˜ − x )u(x ˜u(x) ¯ ¯ ) dV , = G(x (16) ˜ We adopt a Gaussian filter for G(x) and use a filter width of either 2∆ or 4∆, depending on the choice of grid resolution. The SGS heat flux is approximated in the similarity method using (Germano, 1986): Iisim ≈ Csim (u ¯ i T¯ − u˜¯ i T˜¯ ),
(17)
where the coefficient Csim is used to adjust the amplitude of the modeled SGS heat flux. Similar representations are used for the other SGS terms. Explicit filtering under the FEM framework is complicated because each element is connected irregularly to the neighboring elements and physical values are integrated over elements to obtain the filtered values. The ˜ explicit application of G(x) can be avoided by recog˜ ¯ nizing that u(x) depends primarily on the resolved field ¯ u(x) in the vicinity of x. Consequently, we can replace ¯ u(x) in Eq. (16) with u(x) ¯ + ∂u¯ i /∂xj (xj − xj ). Similar substitutions are used for T¯ (x ). Substituting these truncated Taylor expansions into Eq. (17) and carrying out ˜ the filtering with G(x) yields: ∂u¯ i ∂T¯ , ∂xk ∂xk
(12)
Iisim ≈ Csim γ
(13)
˜ where γ is the one-dimensional second moment of G: ˜ γ = x2 G(x) dx dy dz. (19)
the SGS Reynolds stress: Tij = ui uj − u¯ i u¯ j ,
2.2. Nonlinear gradient model
(18)
the SGS Maxwell stress: ¯ iB ¯ j, Mij = Bi Bj − B
(14)
Tijsim ≈ Csim γ
and the SGS induction term: ¯ k ), ai = eijk (uj Bk − u¯ j B
A similar treatment of the other SGS terms gives
(15)
where eijk is the permutation symbol. In order to solve for the filtered fields, we require models for each of these SGS terms. We choose the nonlinear gradient model to approximate equations (12)–(15). A description of the nonlinear gradient model is given in the next section.
∂u¯ i ∂u¯ j , ∂xk ∂xk
(20)
¯j ¯ i ∂B ∂B , ∂xk ∂xk
(21)
Mijsim ≈ Csim γ and
aisim ≈ Csim γeijk
¯k ∂u¯ j ∂B . ∂xl ∂xl
(22)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
111
Fig. 1. Sketch of the relationship between the physical coordinates and the generalized coordinates. Each element in generalized coordinates is a cube with coordinate values ranging from −1 to 1.
There is no need to use the same coefficient Csim for each SGS term, although in this study we set Csim = 1 in Eqs. (18) and (20)–(22). Fig. 2. The grid used the simulation has 64 nodes in each direction.
2.3. Implementation of the nonlinear gradient model in the FEM We implement the nonlinear gradient model in a plane-layer dynamo using the finite-element method (Matsui and Okuda, 2004). An advantage of this general approach is that we can easily adapt the plane-layer model to a spherical shell model. However, one modification is required to implement the nonlinear gradient model on an unstructured grid. In our calculations, we set nodes at the Chebyshev collocation point in the z-direction, and apply the Gaussian filter in the generalized coordinates of the finite-element method (see Fig. 1). Each element in the physical coordinates (x, y, z) is converted to a cubic element in generalized coordinates (ξ, η, ζ), and physical values are interpolated by the linear functions in each element. The derivation of the nonlinear gradient model is now carried out in the generalized coordinate as ¯ ¯f (ξ) = f¯ (ξ 0 ) + ∂f (ξi − ξi0 ) (23) ∂ξi ξ=ξ0 and the corresponding expressions for SGS terms become ∂u¯ i ∂T¯ Ii ≈ Csim γξ , (24) ∂ξk ∂ξk Tij ≈ Csim γξ
∂u¯ i ∂u¯ j , ∂ξk ∂ξk
Mij ≈ Csim γξ
¯ i ∂B ¯j ∂B , ∂ξk ∂ξk
(25) (26)
and ai ≈ Csim γξ eijk
¯k ∂u¯ j ∂B , ∂ξl ∂ξl
(27)
where the coefficient γξ is given by ˜ γξ = ξ 2 G(ξ) dξ dη dζ.
(28)
Comparisons between the predictions of the nonlinear gradient and the similarity models are described in Section 4.1. 2.4. Geometry and boundary conditions We consider a plane layer dynamo to test the SGS model and to investigate the effects of boundaries. As shown in Fig. 2, we use a cubic domain which has unit length in each direction. The top and bottom surfaces are defined by z = ±0.5, and periodic boundary conditions are applied in the horizontal directions. The boundary conditions for u¯ and T¯ on the top and bottom surfaces are u¯ = 0 on z = ±0.5,
(29)
T¯ = 1.0 on z = −0.5,
(30)
and T¯ = 0 on z = 0.5.
(31)
The boundary conditions for the magnetic field in the DNS are prescribed by setting the normal Poynting flux ¯ × B) ¯ · nˆ = 0). Under this condition, only to zero (e.g. (E the normal magnetic field can diffuse out of the domain. The boundary conditions for the vector potential and magnetic field are given by ¯x ¯y ∂A ∂A ¯ z = 0 on z = ±0.5 = =A ∂z ∂z
(32)
112
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
solved, we solve Eqs. (3)–(7) without the SGS terms −∇Tij , Pm−1 E−1 ∇Mij , −∇ · I, and a. The dimensionless numbers in this simulation are chosen to coincide with values adopted in the plane-layer model of Jones and Roberts (2000). We set Pr = 1.0,
(34)
E = 4.0 × 10−4 ,
(35)
Ra = 1.0 × 10 ≈ 4.0Rac ,
(36)
3
Fig. 3. Time evolution of the spatially averaged kinetic and magnetic energies. The kinetic energy is plotted by dotted line, and the magnetic energy is plotted by the solid line.
and ¯x = B ¯ y = 0 on z = ±0.5. B
(33)
3. Results of a direct numerical simulation A fully resolved numerical simulation is evaluated using a finite element mesh with 64 nodes in each direction. Because the fields are well re-
and Pm = 5.0,
(37)
where Rac is the critical Rayleigh number. A large magnetic Prandtl number ensures that the magnetic field is strong. However, this choice means that the magnetic diffusion time is five times longer the viscous diffusion time. The initial condition for velocity and temperature are taken from a simulation of the thermal convection with the same dimensionless numbers but no magnetic field. The initial value for the magnetic field is specified by weak x-component for the vector potential (e.g. ¯ x = 0.001 sin(πz)). A
Fig. 4. Isosurface of z-component of vorticity (upper panels), and intensity of the magnetic field (lower panels). The results for t = 0.6 are given in left panel, and the results for t = 1.9 are given in right panel. We use isosurfaces of ωz = ±3100 for t = 0.6, and ωz = ±2100 for t = 1.9. These amplitudes of vorticity correspond to the 40% of the range of the negative values for each snapshot. The isosurface of B = 13 for t = 0.6 and B = 20 for t = 1.9 is illustrated by red surface. The isosurface of the vorticity is also shown in the low panels with low opacity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
The time evolution of the volume-averaged kinetic energy Ekin and magnetic energy Emag is shown in Fig. 3. The magnetic energy becomes large at t = 1.5, and this intense field causes a small reduction in the kinetic energy. Changes in the fields can be seen in Fig. 4 by comparing snapshots for t = 0.6 and 1.9. (The snapshot at t = 1.9 is chosen to allow sufficient time for the solution to adjust after the transition to the strong field regime.) The top panel of Fig. 4 shows isosurfaces of vorticity, which tend to align with the rotation axis at t = 0.6. This feature becomes less obvious at t = 1.9. Isosurfaces of the magnetic field amplitude are shown in the bottom panel of Fig. 4. Several aspects of the magnetic field are similar in the two snapshots. For example, intense sheet-like structures in the field are generated between intense vortices in the both snapshots. On the other hand, some isosurfaces of the magnetic field can be found inside positive vortices near the boundaries at t = 0.6. This is not observed at t = 1.9. The magnetic field inside the vortices at t = 0.6 is generated by convergence of fluid near the boundary. The magnetic pressure is small at this time so there is little resistance to the convergence of the field lines by the flow. To investigate the length scale of the field patterns, we consider spectra of the Ekin and Emag defined as following: 2 1 2π(kx x+ky )i Ekin (kx , ky ) = dx dy dz , (38) u¯ i e 2 2 1 2π(kx x+ky )i ¯ Emag (kx , ky ) = dx dy dz . (39) Bi e 2 An interesting feature of the solution is that the (kx , ky ) = (0, 0) component of the magnetic field becomes dominant once the magnetic field is strong. Figs. 5 and 6 show the large scale components of the kinetic and magnetic energies. Furthermore, phase angle between the average horizontal magnetic field
and x-axis θ, which is de ¯ x dV / ¯ x2 + B ¯ y2 dV and sin θ = fined by cos θ = B B
¯ y dV / ¯ x2 + B ¯ y2 dV is plotted as a function of B B time. The average horizontal magnetic field rotates about z-axis. The dominant component of the kinetic energy also changes periodically with the phase angle of the horizontal magnetic field. This result indicates that the convection tends to have a larger length scale in the direction of the dominant horizontal magnetic field. Similar behavior was previously observed in the calculations of Jones and Roberts (2000) using comparable dimensionless numbers and the magnetogeostrophic approximation. In this study we retain the full inertia and time step the fluid velocity.
113
Fig. 5. Time evolution of the magnetic (top) and kinetic (middle) energies with (kx , ky ) = (0,1), (1,0), and (1,1) modes and (bottom) the phase angle of the average horizontal magnetic field for 1.8 < t < 2.5.
The behavior of the large-scale components, prior to the increase in magnetic energy, is shown in Fig. 5. It is noted that the averaged horizontal magnetic field also rotates with a period that is slightly longer than in the strong magnetic field case. However, (kx , ky ) = (1, 1) component is the dominant component of the kinetic energy when the magnetic field is weak (see Fig. 6). The (kx , ky ) = (0, 0) component of the magnetic energy is much smaller in Fig. 6, as compared with that in Fig. 5. To investigate the change of the spectra of the energies, the one-dimensional spectra of the kinetic energy are obtained by summing over one of the horizontal wave numbers (e.g. Ekin,x (kx ) = ky Ekin (kx , ky ) and Ekin,y (ky ) = kx Ekin (kx , ky )). One-dimensional spectra of magnetic energies Emag,x (kx ) and Emag,y (ky ) are obtained in the same way. The results for t = 0.6 and 1.9 are plotted in Fig. 7. Changes in the spectra reflect changes in the structure of the fields between the two
114
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
Fig. 7. One-dimensional power spectra of the kinetic energy and the magnetic energy in the resolved DNS.
Fig. 6. Same as Fig. 5 but from the time interval 0.4 < t < 1.1.
snapshots. The large-scale velocities become weaker at t = 1.9, whereas the intermediate scales (10 < k < 20) become slightly stronger. Much of the change in the magnetic field at t = 1.9 occurs at the largest scales. Differences in the kx and ky spectra provide a measure of the anisotropy in the fields. These difference are larger at t = 1.9 when the magnetic field is strong. The relatively rapid decay of the kx spectra (compared with ky ) is due to the orientation of the dominant magnetic field, which is nearly aligned with the x-axis at t = 1.9. 4. Assessing the SGS models
the predictions of the similarity model. The numerical simulation from the previous section is well resolved because the energies at the truncation level are 10−4 times smaller than the largest component of the field (see Fig. 7). We use the solution on a 643 grid to obtain the fields on a coarser 323 grid using a cut-off filter in the spectral domain. We predict the SGS terms in one of three ways: (1) A direct estimation of the SGS terms (Eqs. (12)– (15)) by explicitly filtering the resolved solution on the 643 grid. We use a Gaussian filter with a width of 4∆f , where ∆f refers to the grid spacing on the 643 grid. (2) The similarity model is evaluated by explicitly filtering the fields on the 323 grid. We use a Gaussian filter with a width of 2∆c , where ∆c refers to the grid spacing on the 323 grid. (3) The nonlinear gradient model is evaluated using the fields on the 323 grid.
4.1. Estimates for the SGS terms In this subsection, we evaluate the SGS terms using a snapshot from the direct numerical simulation at t = 1.9. We compare the predicted SGS terms with an estimate obtained from the fully resolved calculation, and with
Although the difference in resolution between the coarse and fine grid is small, it is sufficient for permit an initial test of the SGS models. Successful models at the current resolution do not guarantee good performance in subsequent calculations when the resolution of the
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
115
Fig. 8. Intensity of divergence of the SGS heat flux (upper panels) and z-component of divergence of the SGS momentum flux on the mid-depth of the simulation domain (e.g. z = 0). The result obtained by direct estimation on the resolved grid is given in the left, that obtained by the similarity model is plotted in the middle, and the result of the gradient model is plotted in the right panel.
coarse and fine grids are substantially different. On the other hand, failure at the current resolution would almost certainly ensure failure at more disparate resolutions. The results for the divergence of the SGS heat flux −∇ · I and the z-component of the other SGS terms −∇T , E−1 Pm−1 ∇M, and a at the mid-depth of the domain are shown in Figs. 8 and 9. As seen in Figs. 8 and 9, general features of these three predictions are similar to each other. These comparisons indicate that the nonlinear gradient model is nearly equivalent to the similarity model, and that both models provide a good description of the direct estimation in the mid-plane of the simulation domain. To compare these predictions in more detail, we plot the spectra of the SGS terms with respect to the ydirection in Fig. 10. We can find many common characteristics among all of the SGS terms obtained by the nonlinear gradient model. For example, the SGS terms of the gradient model have larger small-scale components than the predictions of the similarity model. Because of these differences, the small-scale components of the nonlinear gradient model are in better agreement with the small-scale components of the direct estimation.
Only the SGS Lorentz term deviates substantially from the direct estimation. The SGS Lorentz term predicted by the similarity model is even worse. These discrepancies could be reduced by changing the magnitude Csim for the SGS Lorentz term. In general, the magnitude of each SGS term is approximately 10% of the resolved nonlinear terms in the governing equations. One exception occurs in the SGS heat flux, where the SGS term is approximately 4% of the resolved heat flux. The reason of this difference is that the spectrum of the resolved heat flux has a steep slope compared with the spectra of the other resolved terms. We also note that SGS Lorentz force has much larger amplitude than the Reynolds stress. The SGS terms predicted around boundaries are different from those predicted in the interior. The intensity of the divergence of the SGS heat flux and z-component of the SGS Lorentz term on a surface near the lower boundary (z = −0.4975) at t = 1.9 is shown in Fig. 11. The pattern of z-component of SGS Lorentz term is similar to the direct estimation, but the intensity of the SGS terms are smaller than the direct estimates. The divergence of heat flux from the nonlinear gradient model
116
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
Fig. 9. Intensity of z-component of SGS Lorentz term (upper panels) and z-component of SGS induction term (lower panels) term on the mid-depth of the simulation domain (e.g. z = 0). The result obtained by direct estimation on the resolved grid is given in the left, that obtained by the similarity model is plotted in the middle, and the result of the gradient model is plotted in the right panel.
differs in pattern and amplitude from the results of both the direct estimation and similarity model. These differences are most evident in the top-right corner of the plots, where intense negative values are predicted by the nonlinear gradient model, whereas intense positive values are obtained in the direct estimation and similarity model. These results suggest that the SGS terms in the nonlinear gradient model require modification around boundaries. 4.2. Large-eddy simulation using the SGS model As described in the previous section, the SGS terms can be reliably approximated using the nonlinear gradient model in the region away from the boundary layer. To test the utility of this model, we perform a large-eddy simulations with a 323 grid from t = 1.9 to 2.5, using the snapshot at t = 1.9 as an initial value. We compare results obtained from three simulations: (1) A direct simulation with 643 nodes (resolved DNS). (2) A direct simulation with 323 nodes (unresolved DNS).
(3) A simulation using the nonlinear gradient model with 323 nodes (LES). We plot time evolution of the kinetic and magnetic energy averaged over simulation domain for the three cases in Fig. 12. A comparison of the magnetic energy in Fig. 12 shows that the resolved DNS and LES are quite similar prior to t = 2.05. Both of these solutions have large increases in the magnetic energy at t = 2.05, whereas the increase in the unresolved DNS is much smaller. However, the behavior of the LES deviates from the resolved DNS case after t = 2.1. Because of the large time variation in the kinetic and magnetic energies, we compare the time average of the energies from t = 2.2 to 2.5 to avoid the effects of the initial conditions (see Table 1). The average kinetic and magnetic energies in the LES is intermediate between the resolved DNS case and the unresolved DNS case. To compare the convection and magnetic field pattern, we choose a cross section at the mid-plane of the domain after taking 20,000 time steps from the initial value. The corresponding solution time is t = 1.95, which is prior to the divergence of the magnetic energy at t ∼ 2.05. As
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
117
Fig. 10. One-dimensional spectra of SGS terms in the mid-depth of the domain (z = 0.0) at t = 1.9. The prediction of the nonlinear gradient model is plotted by circles with the thick line, while the result of the similarity model is plotted by triangles. The direct estimate on the 643 grid is plotted by ‘+’, and the corresponding resolved terms are plotted by cross with dotted line for reference.
Fig. 11. Intensity of divergence of the SGS heat flux (upper panels) and z-component of the SGS Lorentz term (lower panels) on a surface near the lower boundary (e.g. z = −0.4975). The direct estimation is given in the left, while the result of the similarity model is plotted in the middle. The result of the gradient model is plotted in the right panels.
118
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123 Table 1 Average kinetic and magnetic energies for 2.2 < t < 2.5
Fig. 12. Time evolution of the kinetic energy (top panel) and the magnetic energy (bottom panel) for the resolved DNS, unresolved DNS, and the LES using the nonlinear gradient model.
Run
Resolved DNS
Unresolved DNS
With SGS model
Ekin Emag
7500.3 42922.2
8070.1 49139.3
7727.0 43382.4
seen in Fig. 13, we can find similarities in the pattern of the fields, although the details differ. In general, we can see that the unresolved DNS case has more power in small scale fields. To quantitatively investigate the similarity of the field patterns, we plot the time evolution of the cross correla¯ z in the mid-plane with the resolved DNS tion of the B results in Fig. 14. The correlation coefficient decreases rapidly to 0 by t = 1.95 in both cases, and it is difficult to distinguish between the LES result and unresolved DNS result on the basis of this diagnostic (see Fig. 14). On the other hand, the phase angle between the average horizontal magnetic field and the x-axis is almost same for each model at t = 1.95 (see Fig. 15). These results suggests that the small scale components deviate rapidly because the correlation is strongly influenced by the small scale fields. However, the phase angle also deviates after t = 2.05 (see Fig. 15). The phase angle of the
Fig. 13. Intensity of the z-component of vorticity (upper panels), and magnetic field (lower panels) in the mid-plane of the domain (z = 0.0) at t = 1.95. The result obtained by the resolved DNS is given in the left, the result obtained by the unresolved DNS is plotted in the middle, and the result by the LES is plotted in the right panel.
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
119
Fig. 14. Time evolution of the correlation coefficient for the zcomponent of the magnetic field in the mid-plane. The correlation between the resolved and unresolved DNS is plotted by dotted line, and the correlation between the resolved DNS and LES is plotted by solid line.
LES rotates more slowly than that in the resolved DNS case, whereas the phase angle in unresolved DNS case changes faster than that in the resolved DNS result. To investigate the effects of the SGS terms on the resolved fields, we plot a power spectrum of the kinetic and magnetic energies with respect to the y-direction for t = 2.4 in Fig. 16. The spectrum of the kinetic energy in the LES is much closer to the kinetic energy spectrum in resolved direct simulation. The unresolved DNS substantially overestimates the kinetic energy at the small scales. This result suggests that the SGS terms provide adequate dissipation of the kinetic energy at the small scales. The spectra of magnetic energy have similar characteristics, but the small scale components of the magnetic energy in the LES become much smaller than those in the resolved DNS case. These characteristics also can be seen the snapshot of the fields at t = 1.95 (see Fig. 13). Recalling the average kinetic and magnetic energies in Table 1, the difference of the averaged energy is due mainly to the difference in the large scale fields. The relative importance of the SGS terms in the momentum equation can be inferred from the power spectra
Fig. 15. Time evolution of the phase angle of the average horizontal magnetic field for the resolved DNS (thick solid line), unresolved DNS (dashed line), and the LES using the nonlinear gradient model (thin solid line).
Fig. 16. One-dimensional spectra of the kinetic and magnetic energies at t = 2.4.
of the SGS terms (see Fig. 10). The SGS Lorentz term is approximately 10 times larger than the Reynolds stress. In fact, the SGS Lorentz term is comparable to the resolved inertial term. The magnitudes of the SGS Lorentz and Reynolds stress suggest that the Reynolds stress has a smaller role in the present simulation, and that the difference in the kinetic energy spectrum between the LES and unresolved DNS cases is due to the SGS Lorentz term. Energy balances are also important for understanding the role of the SGS terms. The conservation of the resolved kinetic energy is obtained by integrating the product of u¯ and Eq. (3) over the fluid volume. The resolved magnetic energy can similarly be described by ¯ and Eq. (5) over the fluid integrating the product of B volume. The result is ∂ 1 2 u¯ dV ∂t 2 1 2 =− u¯ u¯ + u¯ χ¯ + ω¯ × u¯ · dS 2 ¯ + −ω¯ 2 + RaE−1 u¯ z (T¯ − T0 ) + Pm−1 E−1 u· ¯ − u¯ · ∇Tij + Pm−1 E−1 u¯ · ∇Mij dV, (J¯ × B) (40)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
120
Fig. 17. Time evolution of energy flux for the kinetic (upper panel) and magnetic energies (lower panel) for the LES. It is noted that all of the SGS contributions are negative. The work done by the Reynolds stress is not shown in the upper panel because it is not distinguishable from zero in this plot.
∂ ∂t
1 ¯2 B dV 2
=−
¯ ×B ¯ dS E
[−Pm−1 J¯ 2
¯ − J¯ · a] dV, − J¯ · (u¯ × B)
(41)
Fig. 18. Schematic illustration of the role of the SGS induction and Lorentz terms in the kinetic and magnetic energies. The work done by the resolved Lorentz force transfers kinetic energy from the resolved velocity into the resolved magnetic energy. Similarly, the SGS Lorentz term draws kinetic energy from the resolved velocity into magnetic energy of the unresolved field. Generation of magnetic energy by the resolved fields is positive, whereas the SGS induction term dissipates magnetic energy.
dicating a sink of resolved magnetic energy. By contrast the resolved induction term is positive, indicating a constructive contribution to the magnetic energy. We summarize the energy transfer between resolved and unresolved component by the SGS terms in Fig. 18. The SGS Lorentz term and SGS induction term transfer energy from the resolved components of the kinetic and magnetic energies into the unresolved components.
where ¯ = Pm−1 J¯ − (u¯ × B) ¯ − a. E
(42)
The boundary conditions (Eqs. (30)–(33)) ensure that the surface integrals in Eqs. (40) and (41) vanish in the present model. The importance of the various terms in Eqs. (40) and (41) is shown in Fig. 17. The SGS Lorentz force and induction term have significant role in the energy balance because they are large compared with the other SGS terms. As seen in Fig. 17, the work done by the resolved and SGS Lorentz force has same sign through the simulation. On the other hand, the SGS induction term is negative throughout the simulation in-
5. Discussion The results described in the previous sections indicate that the SGS model transfers energy from large to small scales. The SGS Lorentz term is primarily responsible for dissipating kinetic energy in the resolved scales, whereas the SGS induction term dissipates large-scale magnetic energy. Comparison of the magnetic energy spectra in Fig. 16 suggests that the SGS induction term draws too much energy from the resolved components. The fact that the nonlinear gradient model reproduces the spectral amplitude of the direct estimate in Fig. 10 suggests
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
121
Fig. 19. Intensity of the y-component (upper panels) and z-component (lower panels) of the SGS heat flux on a surface near the lower boundary (e.g. z = −0.4975). The direct estimation is given in the left, while the result of the similarity model is plotted in the middle. The result of the gradient model is plotted in the right panels.
that the phase of those two terms is slightly different. One way to reduce the dissipation caused by the SGS induction term is to reduce the model constant Csim in Eq. (27), although the spectral amplitudes of the SGS induction term would then be small. Similar problems arise with the SGS Lorentz term. The amplitude of the SGS Lorentz term is much smaller than that predicted by the direct estimation (see Fig. 10), yet the dissipation of the kinetic energy is about right. This inconsistency may be a consequence of comparing spectra at a specific snapshot in time. The spectra of the SGS terms relate to the time derivative of the energies, rather than the energies directly. It would be preferable to base comparisons on time-averaged spectra. Another problem with the current implementation of the nonlinear gradient model occurs near the boundaries. We previously noted in Fig. 11 that the predicted pattern of −∇ · I near the boundary differs from that of the direct estimation. In order to identify the origin of this discrepancy, we plot the components Iy and Iz near the boundary at the same time step. We can see that these components are in much better agreement with the SGS heat flux from the direct estimation, although the
amplitudes are slightly different. The nonlinear gradient model appears to overestimate Iy and underestimate Iz . If we repeat this comparison at different distances from the boundaries, we find that the discrepancy in the amplitude of the heat flux components varies systematically with distance. Specifically, the nonlinear gradient model increasingly overestimates the direct estimate as we approach the boundary. This suggests that the model coefficients Csim should decrease toward the boundary. Germano et al. (1991) devised a scheme to obtain the model coefficient Csim by explicitly filtering the resolved fields. This scheme was originally developed to estimate the coefficient for the Smagorinsky model. The success of the so-called dynamic model was demonstrated by Piomelli (1993), who showed that the coefficient of the Smagorinsky model approaches to 0 toward the boundary. The dynamic model is quite general and can be applied to other SGS models, including the nonlinear gradient model. Thus, the potential exists to evaluate the coefficients of the nonlinear gradient model using a dynamic scheme. If this approach can be used to approximate the constant Csim (z) that we obtained from a fit to the direct estimation of the SGS terms, then we
122
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
Fig. 20. Comparison of the SGS heat flux when calculated using −∇ · (uT ) + ∇ · (u¯ T¯ ) (left panel) and −[∇ · (uT − u¯ T¯ )] (middle panel). Both terms are evaluated on the 643 grid near the lower boundary (i.e. z = −0.4975). The result of the nonlinear gradient model on the 323 grid is plotted in the right panel for comparison.
can expect good performance in large-eddy simulations of convection-driven dynamos. Another source of difficulty near the boundaries may be caused by the use of a non-uniform grid. We have seen in Fig. 19 that the pattern of Iy and Iz of the gradient model around the boundary is quite similar to that of the direct estimation and the similarity model. However, the pattern of −∇ · I by the gradient model is quite different from that obtained by the direct estimation and similarity model (see Fig. 11). Ghosal and Moin (1995) pointed out that the filtering operation does not, in general, commute with the operation of the differentiation if we use a filter with a variable filter width. Commutation of the operations is customarily assumed in deriving the governing equations for the large-scale fields. For example, in the heat equation for the large scale temperature field, we assume that ∂ ∂ (ui T ) = (ui T ). ∂xi ∂xi
(43)
In the present study, the size of element changes with respect to the horizontal position (see Fig. 2) and filter width is changed by the size of the element in the physical space. To asses the error caused by interchanging the order of differentiation and filtering, we compare two representations of the SGS heat flow. One representation is given by −∇ · (uT ) + ∇ · (u¯ T¯ ), which was previously used to define the direct estimate in Fig. 11. The other representation is given by −[∇ · (uT − u¯ T¯ )]. Both terms are evaluate near the boundary, where the grid spacing is most uneven. Differences between the two terms provide a measure of the commutation error Fig. 20. We note that the pattern of −[∇ · (uT − u¯ T¯ )] is similar to −∇ · I sim , as obtained using the nonlinear gradient model. Consequently the error in the SGS model near the boundaries
may be a result of the commutation error. Vasilyev et al. proposed an approach to construct a commutative filter for the unstructured grid. Construction of a commutative filter will be important for improving the SGS model especially near the boundaries. 6. Conclusion We implemented a sub-grid scale model (SGS) in a dynamo simulation code using the finite-element method. We applied the nonlinear gradient model which is a version of the similarity SGS model to the four nonlinear terms in the MHD equations, e.g. the heat and momentum fluxes, Lorentz term, and the induction terms. First, we test the predictions of the SGS model using snapshots from a direct numerical simulation (DNS). The result showed that the spatial pattern of the SGS terms in the nonlinear gradient model is similar to that obtained by the direct estimation from a resolved DNS at the mid-plane. The only exception occurs with SGS Lorentz term, where the nonlinear gradient model underestimates the direct estimation. Simulations using this SGS model are performed in a rotating plane layer model, and we compared the results with the output a direct numerical simulation using either the same or finer spatial resolution. All of those simulations have a large time variation in the kinetic and magnetic energies. However, the time averaged kinetic and magnetic energies energy obtained by LES are intermediate between the resolved DNS case and the unresolved DNS case. Comparing field patterns, the patterns of convection and magnetic field in the LES are more like those in the resolved DNS case than those in the unresolved DNS because the unresolved results have strong small-scale fields.
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
We also investigate the energetics of the convection and magnetic field in the simulation. The work done by the SGS Lorentz term and the magnetic energy generated by the SGS induced terms are both negative, indicating that energy is transferred form the resolved components to the unresolved components. While the simulations suggests that the SGS models works adequately, we also find some discrepancy between the LES and resolved DNS results. The spectrum of the kinetic energy in the LES is much closer to that in the resolved DNS (as compared with the unresolved DNS). On the other hand, small-scale components in the magnetic energy of the LES are smaller than those in the both resolved and unresolved DNS, implying that the SGS model is overly dissipative. The largest discrepancies in the SGS models occur near the boundaries. To improve the present model, we need to modify the filter to reduce the commutation error caused by interchanging the order of filtering and spatial differentiation. The dynamic model also has the potential to improve the prediction of SGS terms around the boundaries by applying the model coefficients as a function of the vertical position z. Acknowledgements The work is supported by the National Science Foundation. The authors would like to the National institute for Polar Research and to Center for Planning and Information systems in Japan Aerospace Exploration Agency for performing the present simulations. References Braginsky, S.I., Meytlis, V.P., 1990. Local turbulence in the earth’s core. Geophys. Astrophys. Fluid Dyn. 55, 71–87. 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. Buffett, B.A., 2003. A comparison of subgrid-scale models for largeeddy simulations of convection in the earth’s core. Geophys. J. Int. 153, 753–765. Christensen, U.R., Olson, P., Glatzmaier, G.A., 1999. Numerical modelling of the geodynamo: a systematic parameter study. Geophys. J. Int. 138, 393–409. Germano, M., 1986. A proposal for a redefinition of the turbulent stresses in the filtered navier stokes equations. Phys. Fluids 29, 2323–2324. Germano, M., Piomelli, U., Moin, P., Cabot, W., 1991. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids 3, 1760–1765. Ghosal, S., Moin, P., 1995. The basic equations for the large eddy simulation of turbulent flows in complex geometry. J. Comput. Phys. 118, 24–37.
123
Glatzmaier, G.A., Roberts, P.H., 1995. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle. Phys. Earth Planet. Interiors 91, 63–75. Glatzmaier, G.A., Roberts, P.H., 1997. Simulating the geodynamo. Contemp. Phys. 38, 269–288. Jones, C.A., Roberts, P.H., 2000. Convection-driven dynamos in a rotating plane layer. J. Fluid Mech. 404, 311–343. Kageyama, A., Sato, T., the Complexity Simulation Group, 1995. Computer simulation of a magnetohydrodynamic dynamo. Phys. Plasmas 2, 1421–1431. Kuang, W., Bloxham, J., 1997. An earth-like numerical dynamo model. Nature 389, 371–374. Kuang, W., Bloxham, J., 1999. Numerical modeling of magnetohydrodynamic convection in a rapidly rotating spherical shell: weak and strong field dynamo action. J. Comput. Phys. 153, 51–81. Leonard, A., 1974. Energy cascade in large-eddy simulations of turbulent fluid flows. Adv. Geophys. 18, 237–248. Lesieur, M., Metais, O., 1996. New trends in large-eddy simulations of turbulence. Annu. Rev. Fluid Mech. 28, 45–82. Matsui, H., Okuda, H., 2004. Development of a simulation code for mhd dynamo processes using the geofem platform. Int. J. Comput. Fluid Dyn. 18 (4), 323–332. Matsushima, M., 2001. Expression of turbulent heat flux in the earth’s core in terms of a second moment closure model. Phys. Earth Planet. Interiors 104, 137–148. Matsushima, M., Nakajima, T., Roberts, P.H., 1999. The anisotoropy of local turbulence in the earth’s core. Earth Planets Space 51, 277–286. Meneveau, C., Katz, J., 2000. Scale-invariance and turbulence models for large-eddy simulation. Annu. Rev. Fluid Mech. 32, 1–32. Moffatt, H.K., Proctor, M.R.E., 1982. The role of the helicity spectrum function in turbulent dynamo theory. Geophys. Astrophys. Fluid Dyn. 21, 265–283. Parker, E.N., 1955. Hydromagnetic dynamo models. Astrophys. J. 122, 293–314. Phillips, C.G., Ivers, D.J., 2001. Special interactions of rapidly-rotating anisotropic turbulent viscous and thermal diffusion in the earth’s core. Phys. Earth Planet. Interiors 128, 93–107. Phillips, C.G., Ivers, D.J., 2003. Strong field anisotropic diffusion models for the earth’s core. Phys. Earth Planet. Interiors 140, 13–28. Piomelli, U., 1993. High reynolds number calculations using the dynamic subgrid-scale stress model. Phys. Fluids A 5, 1484– 1490. R¨adler, K.-H., Wiedemann, E., Brandenburg, A., Meinel, R., Tuominen, I., 1990. Nonlinear mean-field dynamo models: stability and evolution of three-dimensional magnetic field configurations. Astron. Astrophys. 239, 413–423. Sakuraba, A., Kono, M., 1999. Effects of the inner core on the numerical simulation of the magnetohydrodynamic dynamo. Phys. Earth Planet. Interiors 111, 105–121. Shimizu, H., Loper, D.E., 1997. Time and length scales of buoyancydriven flow structures in a rotating hydromagnetic fluid. Phys. Earth Planet. Interiors 104, 307–329. Winckelmans, G.S., Wray, A.S., Vasilyev, O.V., Jeanmart, H., 2001. Explicit-filtering large-eddy simulation using the tensor-diffusivity model supplemented by a dynamic smagorinsky term. Phys. Fluids 13, 1385–1403. Zhang, K.K., Jones, C.A., 1997. The effect of hyperviscosity on geodynamo models. Geophys. Res. Lett. 24 (22), 2869–2872.
Sub-grid scale model for convection-driven dynamos in a rotating plane layer Hiroaki Matsui∗ , Bruce A. Buffett Department of Geophysical Sciences, The University of Chicago, 5734 S Ellis Avenue, Chicago, IL 60637, USA Received 2 November 2004; accepted 3 March 2005
Abstract We develop a sub-grid scale (SGS) model for convection-driven dynamos, based on the nonlinear gradient model of Leonard [Leonard, A., 1974. Energy cascade in large-eddy simulations of turbulent fluid flows. Adv. Geophys. 18, 237–248]. The predictions of the SGS model are tested using snapshots from a direct numerical simulation (DNS). Good agreement is obtained throughout most of the fluid, but discrepancies occur near the boundaries. We also implement the SGS model in a large-eddy simulation (LES) to assess the temporal behavior of the SGS model. Both the DNS and LES have large time variation in the kinetic and magnetic energies, so comparisons are made using time averages. The time-averaged energies in the LES are slightly larger than those in the resolved DNS, but not as high as those obtained in an unresolved DNS with the same grid resolutions as the LES. The source of dissipation in the LES is revealed by examining the energy balances. The work done by the SGS Lorentz term and the magnetic energy generated by the SGS induction term are both negative, indicating that these SGS terms transfer energy from the resolved scales into the unresolved scales. © 2005 Elsevier B.V. All rights reserved. Keywords: Geodynamo; Magnetohydrodynamics; Sub-grid scale model
1. Introduction In the last 10 years, geodynamo simulations have begun to reveal the basic processes responsible for generating the geomagnetic field. Following the initial studies of Glatzmaier and Roberts (1995) and Kageyama and Sato (1995), a number of additional geodynamo models have been developed, which reproduce many characteristics of the geomagnetic field (e.g. Kuang and Bloxham, 1997, 1999; Christensen et al., 1999; Sakuraba and Kono, 1999). However, the fluid motions in the outer core are expected to have a broad range of length scales, ranging from the dissipative scale in the boundary layer (less than L ∼ 10 m) to the geometry of the outer core ∗
Corresponding author. Tel.: +1 773702 8107; fax: +1 773702 9505. E-mail address: [email protected] (H. Matsui)
0031-9201/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2005.03.019
(L ∼ 106 m). Current numerical simulations cannot use realistic parameters for the Earth’s core because of limited spatial resolution. To account for the influence of unresolved flow, all simulations have so far used turbulent diffusivities (e.g. Glatzmaier and Roberts, 1995; Sakuraba and Kono, 1999). Some of these studies use constant diffusivities, whereas others allow the diffusivities to depend on the length scale of the motion (e.g. hyperdiffusivity). Zhang and Jones (1997) noted that hyperdiffusivity alters the structure of large scale convection, which means that the choice of sub-grid scale models can affect the resolved fields. Small scale convection in the outer core is expected to be strongly anisotropic due to the influences of the Earth’s rotation and magnetic field (Braginsky and Meytlis, 1990; Shimizu and Loper, 1997; Matsushima et al., 1999). However, few models deal with anisotropic
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
diffusivities in numerical calculations. One exception is described by Phillips and Ivers (2001, 2003), who have developed a tensor diffusivity model for use in pseudospectral geodynamo codes. Another example is given by Matsushima (2001), who used magnetoconvection calculation in a small periodic domain of the outer core to devise an anisotropic model for turbulent thermal diffusion. Buffett (2003) also used magnetoconvection calculations to test predictions of SGS heat and momentum fluxes using turbulent viscosity, hyper diffusivity, as well as the Smagorinsky and the similarity models. He concluded that the similarity model successfully reproduces the anisotropy of the SGS heat flux with respect to the rotation axis and direction of the imposed magnetic field. In the present study, we implement a variant of the similarity model in a MHD dynamo model using the finite-element method (FEM). The model is based on the nonlinear gradient model of Leonard (1974), which represents the SGS terms using the gradients of the resolved fields. This approach does not require explicit filtering of the resolved field, as is often the case in SGS modeling (see Lesieur and Metais, 1996; Meneveau and Katz, 2000 for reviews). We apply the nonlinear gradient model to a plane-layer dynamo. Large-eddy simulations (LES) are performed using the nonlinear gradient model and the results are compared with estimates from a direct numerical simulation (DNS) using a finer grid. The goals of the present study are three-fold. First, we attempt to identify and develop practical methods that are capable of reproducing the strong anisotropy expected in the SGS heat and momentum flux. Encouraging results have previously been obtained using the similarity model in magnetoconvection calculations (Buffett, 2003; Matsushima, 2001). However, methods have not yet been tested in dynamo calculations. Second, we need to extend the similarity model to include SGS Lorentz and induction terms. Previous dynamo models in terrestrial (e.g. Glatzmaier and Roberts, 1997) and stellar contexts (e.g. Brun et al., 2004) account for only the eddy diffusivities. The extension of the similarity model to include the SGS Lorentz and induction terms is straightforward, and our calculations show that these additional terms contribute significantly to the evolution of the kinetic and magnetic energy at the resolved scales. In fact, the similarity model can potentially transfer energy into the resolved scales from the unresolved scales. Such a process might emulate the α-effect, which was originally proposed to describe the generation the large scale magnetic field by small scale helical flow (e.g. Parker, 1955; Moffatt and Proctor, 1982; R¨adler et al., 1990). However, the transfer of the energy into the resolved fields may also causes numerical instabilities. Previous stud-
109
ies of homogeneous turbulence have shown that the similarity model does not dissipate enough kinetic energy (e.g. Winckelmans et al., 2001). We test the similarity model in the dynamo problem without the addition of an eddy diffusivity term to address this specific question. We show that the SGS model may actually overestimate the dissipation of magnetic energy in the resolved scales. This and other shortcomings motivate proposals for further improvements to the SGS similarity model. 2. Simulation models 2.1. Sub-grid scale modeling Sub-grid scale modeling often relies on spatial filtering to eliminate field components with length scales smaller than the grid spacing ∆. Each field in the model (e.g. velocity u, . . .) is convolved with a filter function G(x) to define the large scale field (e.g. Leonard, 1974): u¯ = G(x − x )u(x ) dV, (1) where u¯ is called the filtered or resolved velocity. Many types of filter functions are available, and we adopt a Gaussian filter, given by 6 1 6x2 G(x) = exp − 2 . (2) π∆ ∆ The filter is applied to the governing equations to obtain a set of equations for resolved fields. The filtering procedure introduces one new term into the governing equations for each nonlinear term in the original equations (see, for example, Buffett, 2003). We consider the problem of convection-driven motion of the conductive fluid in a plane layer that rotates with the angular velocity = Ωˆz. We use the Boussinesq approximation and take gravity to be a constant vector g = −gˆz. The fluid has a kinetic viscosity ν, thermal diffusivity κ, thermal expansion coefficient α, and magnetic diffusivity η. To obtain the normalized equations, we use the thickness D of the fluid layer as a length scale, and the viscous diffusion time, D2 /ν, as the characteristic time scale. The magnetic field is normalized by B02 = ρµ0 ηΩ, where ρ is the density of the fluid and µ0 is the magnetic permeability. We apply the filter to the basic equations for a conducting fluid with the Boussinesq approximation to obtain ∂u¯ ¯ + ∇(u¯ u) ∂t ¯ + RaE−1 (T¯ − T0 )ˆz = −∇ χ¯ + ∇ 2 u¯ − 2E−1 (ˆz × u) −1 −1 ¯ B)−∇T ¯ + Pm−1 E−1 ∇(B ij +Pm E ∇Mij , (3)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
110
∂T¯ = −∇ · (u¯ T¯ ) + Pr −1 ∇ 2 T¯ − ∇ · I, ∂t
(4)
¯ ∂A ¯ + u¯ × B ¯ + a, = −∇ ϕ¯ + Pm−1 ∇ 2 A ∂t
(5)
¯ = 0, ∇ · u¯ = ∇ · A
(6)
and ¯ ¯ = ∇ × A, B
(7)
where u, χ, T, T0 , B, A, and ϕ are the velocity, modified pressure, temperature perturbation, reference temperature, magnetic field, vector potential of the magnetic field, and arbitrary scalar potential to satisfy the Coulomb gauge, respectively. The bar indicates the filtered fields. The reference temperature T0 is given by T0 = (ztop − z)/D, where z is the vertical coordinate. The Prandtl number Pr, the modified Rayleigh number Ra, the Ekman number E, and the magnetic Prandtl number Pm are defined as follows: ν Pr = , (8) κ αg∆TD , Ων ν E= , ΩD2 Ra =
(9) (10)
and Pm =
ν . η
(11)
Four SGS terms appear in the basic equations (3)–(5), each corresponding to a nonlinear term in the unfiltered equations. These include the SGS heat flux: Ii = ui T − u¯ i T¯ ,
The nonlinear gradient model was original proposed by Leonard (1974). It can now be viewed as a version of the similarity model, which is based on the scale invariance of the turbulence. We describe this model using the SGS heat flux I = uT − u¯ T¯ as an example. The similarity model relies on a second filter with a width that ˜ is greater than the grid size ∆. The second filter G(x) is ¯ applied to the resolved field u(x) to give ˜ − x )u(x ˜u(x) ¯ ¯ ) dV , = G(x (16) ˜ We adopt a Gaussian filter for G(x) and use a filter width of either 2∆ or 4∆, depending on the choice of grid resolution. The SGS heat flux is approximated in the similarity method using (Germano, 1986): Iisim ≈ Csim (u ¯ i T¯ − u˜¯ i T˜¯ ),
(17)
where the coefficient Csim is used to adjust the amplitude of the modeled SGS heat flux. Similar representations are used for the other SGS terms. Explicit filtering under the FEM framework is complicated because each element is connected irregularly to the neighboring elements and physical values are integrated over elements to obtain the filtered values. The ˜ explicit application of G(x) can be avoided by recog˜ ¯ nizing that u(x) depends primarily on the resolved field ¯ u(x) in the vicinity of x. Consequently, we can replace ¯ u(x) in Eq. (16) with u(x) ¯ + ∂u¯ i /∂xj (xj − xj ). Similar substitutions are used for T¯ (x ). Substituting these truncated Taylor expansions into Eq. (17) and carrying out ˜ the filtering with G(x) yields: ∂u¯ i ∂T¯ , ∂xk ∂xk
(12)
Iisim ≈ Csim γ
(13)
˜ where γ is the one-dimensional second moment of G: ˜ γ = x2 G(x) dx dy dz. (19)
the SGS Reynolds stress: Tij = ui uj − u¯ i u¯ j ,
2.2. Nonlinear gradient model
(18)
the SGS Maxwell stress: ¯ iB ¯ j, Mij = Bi Bj − B
(14)
Tijsim ≈ Csim γ
and the SGS induction term: ¯ k ), ai = eijk (uj Bk − u¯ j B
A similar treatment of the other SGS terms gives
(15)
where eijk is the permutation symbol. In order to solve for the filtered fields, we require models for each of these SGS terms. We choose the nonlinear gradient model to approximate equations (12)–(15). A description of the nonlinear gradient model is given in the next section.
∂u¯ i ∂u¯ j , ∂xk ∂xk
(20)
¯j ¯ i ∂B ∂B , ∂xk ∂xk
(21)
Mijsim ≈ Csim γ and
aisim ≈ Csim γeijk
¯k ∂u¯ j ∂B . ∂xl ∂xl
(22)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
111
Fig. 1. Sketch of the relationship between the physical coordinates and the generalized coordinates. Each element in generalized coordinates is a cube with coordinate values ranging from −1 to 1.
There is no need to use the same coefficient Csim for each SGS term, although in this study we set Csim = 1 in Eqs. (18) and (20)–(22). Fig. 2. The grid used the simulation has 64 nodes in each direction.
2.3. Implementation of the nonlinear gradient model in the FEM We implement the nonlinear gradient model in a plane-layer dynamo using the finite-element method (Matsui and Okuda, 2004). An advantage of this general approach is that we can easily adapt the plane-layer model to a spherical shell model. However, one modification is required to implement the nonlinear gradient model on an unstructured grid. In our calculations, we set nodes at the Chebyshev collocation point in the z-direction, and apply the Gaussian filter in the generalized coordinates of the finite-element method (see Fig. 1). Each element in the physical coordinates (x, y, z) is converted to a cubic element in generalized coordinates (ξ, η, ζ), and physical values are interpolated by the linear functions in each element. The derivation of the nonlinear gradient model is now carried out in the generalized coordinate as ¯ ¯f (ξ) = f¯ (ξ 0 ) + ∂f (ξi − ξi0 ) (23) ∂ξi ξ=ξ0 and the corresponding expressions for SGS terms become ∂u¯ i ∂T¯ Ii ≈ Csim γξ , (24) ∂ξk ∂ξk Tij ≈ Csim γξ
∂u¯ i ∂u¯ j , ∂ξk ∂ξk
Mij ≈ Csim γξ
¯ i ∂B ¯j ∂B , ∂ξk ∂ξk
(25) (26)
and ai ≈ Csim γξ eijk
¯k ∂u¯ j ∂B , ∂ξl ∂ξl
(27)
where the coefficient γξ is given by ˜ γξ = ξ 2 G(ξ) dξ dη dζ.
(28)
Comparisons between the predictions of the nonlinear gradient and the similarity models are described in Section 4.1. 2.4. Geometry and boundary conditions We consider a plane layer dynamo to test the SGS model and to investigate the effects of boundaries. As shown in Fig. 2, we use a cubic domain which has unit length in each direction. The top and bottom surfaces are defined by z = ±0.5, and periodic boundary conditions are applied in the horizontal directions. The boundary conditions for u¯ and T¯ on the top and bottom surfaces are u¯ = 0 on z = ±0.5,
(29)
T¯ = 1.0 on z = −0.5,
(30)
and T¯ = 0 on z = 0.5.
(31)
The boundary conditions for the magnetic field in the DNS are prescribed by setting the normal Poynting flux ¯ × B) ¯ · nˆ = 0). Under this condition, only to zero (e.g. (E the normal magnetic field can diffuse out of the domain. The boundary conditions for the vector potential and magnetic field are given by ¯x ¯y ∂A ∂A ¯ z = 0 on z = ±0.5 = =A ∂z ∂z
(32)
112
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
solved, we solve Eqs. (3)–(7) without the SGS terms −∇Tij , Pm−1 E−1 ∇Mij , −∇ · I, and a. The dimensionless numbers in this simulation are chosen to coincide with values adopted in the plane-layer model of Jones and Roberts (2000). We set Pr = 1.0,
(34)
E = 4.0 × 10−4 ,
(35)
Ra = 1.0 × 10 ≈ 4.0Rac ,
(36)
3
Fig. 3. Time evolution of the spatially averaged kinetic and magnetic energies. The kinetic energy is plotted by dotted line, and the magnetic energy is plotted by the solid line.
and ¯x = B ¯ y = 0 on z = ±0.5. B
(33)
3. Results of a direct numerical simulation A fully resolved numerical simulation is evaluated using a finite element mesh with 64 nodes in each direction. Because the fields are well re-
and Pm = 5.0,
(37)
where Rac is the critical Rayleigh number. A large magnetic Prandtl number ensures that the magnetic field is strong. However, this choice means that the magnetic diffusion time is five times longer the viscous diffusion time. The initial condition for velocity and temperature are taken from a simulation of the thermal convection with the same dimensionless numbers but no magnetic field. The initial value for the magnetic field is specified by weak x-component for the vector potential (e.g. ¯ x = 0.001 sin(πz)). A
Fig. 4. Isosurface of z-component of vorticity (upper panels), and intensity of the magnetic field (lower panels). The results for t = 0.6 are given in left panel, and the results for t = 1.9 are given in right panel. We use isosurfaces of ωz = ±3100 for t = 0.6, and ωz = ±2100 for t = 1.9. These amplitudes of vorticity correspond to the 40% of the range of the negative values for each snapshot. The isosurface of B = 13 for t = 0.6 and B = 20 for t = 1.9 is illustrated by red surface. The isosurface of the vorticity is also shown in the low panels with low opacity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
The time evolution of the volume-averaged kinetic energy Ekin and magnetic energy Emag is shown in Fig. 3. The magnetic energy becomes large at t = 1.5, and this intense field causes a small reduction in the kinetic energy. Changes in the fields can be seen in Fig. 4 by comparing snapshots for t = 0.6 and 1.9. (The snapshot at t = 1.9 is chosen to allow sufficient time for the solution to adjust after the transition to the strong field regime.) The top panel of Fig. 4 shows isosurfaces of vorticity, which tend to align with the rotation axis at t = 0.6. This feature becomes less obvious at t = 1.9. Isosurfaces of the magnetic field amplitude are shown in the bottom panel of Fig. 4. Several aspects of the magnetic field are similar in the two snapshots. For example, intense sheet-like structures in the field are generated between intense vortices in the both snapshots. On the other hand, some isosurfaces of the magnetic field can be found inside positive vortices near the boundaries at t = 0.6. This is not observed at t = 1.9. The magnetic field inside the vortices at t = 0.6 is generated by convergence of fluid near the boundary. The magnetic pressure is small at this time so there is little resistance to the convergence of the field lines by the flow. To investigate the length scale of the field patterns, we consider spectra of the Ekin and Emag defined as following: 2 1 2π(kx x+ky )i Ekin (kx , ky ) = dx dy dz , (38) u¯ i e 2 2 1 2π(kx x+ky )i ¯ Emag (kx , ky ) = dx dy dz . (39) Bi e 2 An interesting feature of the solution is that the (kx , ky ) = (0, 0) component of the magnetic field becomes dominant once the magnetic field is strong. Figs. 5 and 6 show the large scale components of the kinetic and magnetic energies. Furthermore, phase angle between the average horizontal magnetic field
and x-axis θ, which is de ¯ x dV / ¯ x2 + B ¯ y2 dV and sin θ = fined by cos θ = B B
¯ y dV / ¯ x2 + B ¯ y2 dV is plotted as a function of B B time. The average horizontal magnetic field rotates about z-axis. The dominant component of the kinetic energy also changes periodically with the phase angle of the horizontal magnetic field. This result indicates that the convection tends to have a larger length scale in the direction of the dominant horizontal magnetic field. Similar behavior was previously observed in the calculations of Jones and Roberts (2000) using comparable dimensionless numbers and the magnetogeostrophic approximation. In this study we retain the full inertia and time step the fluid velocity.
113
Fig. 5. Time evolution of the magnetic (top) and kinetic (middle) energies with (kx , ky ) = (0,1), (1,0), and (1,1) modes and (bottom) the phase angle of the average horizontal magnetic field for 1.8 < t < 2.5.
The behavior of the large-scale components, prior to the increase in magnetic energy, is shown in Fig. 5. It is noted that the averaged horizontal magnetic field also rotates with a period that is slightly longer than in the strong magnetic field case. However, (kx , ky ) = (1, 1) component is the dominant component of the kinetic energy when the magnetic field is weak (see Fig. 6). The (kx , ky ) = (0, 0) component of the magnetic energy is much smaller in Fig. 6, as compared with that in Fig. 5. To investigate the change of the spectra of the energies, the one-dimensional spectra of the kinetic energy are obtained by summing over one of the horizontal wave numbers (e.g. Ekin,x (kx ) = ky Ekin (kx , ky ) and Ekin,y (ky ) = kx Ekin (kx , ky )). One-dimensional spectra of magnetic energies Emag,x (kx ) and Emag,y (ky ) are obtained in the same way. The results for t = 0.6 and 1.9 are plotted in Fig. 7. Changes in the spectra reflect changes in the structure of the fields between the two
114
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
Fig. 7. One-dimensional power spectra of the kinetic energy and the magnetic energy in the resolved DNS.
Fig. 6. Same as Fig. 5 but from the time interval 0.4 < t < 1.1.
snapshots. The large-scale velocities become weaker at t = 1.9, whereas the intermediate scales (10 < k < 20) become slightly stronger. Much of the change in the magnetic field at t = 1.9 occurs at the largest scales. Differences in the kx and ky spectra provide a measure of the anisotropy in the fields. These difference are larger at t = 1.9 when the magnetic field is strong. The relatively rapid decay of the kx spectra (compared with ky ) is due to the orientation of the dominant magnetic field, which is nearly aligned with the x-axis at t = 1.9. 4. Assessing the SGS models
the predictions of the similarity model. The numerical simulation from the previous section is well resolved because the energies at the truncation level are 10−4 times smaller than the largest component of the field (see Fig. 7). We use the solution on a 643 grid to obtain the fields on a coarser 323 grid using a cut-off filter in the spectral domain. We predict the SGS terms in one of three ways: (1) A direct estimation of the SGS terms (Eqs. (12)– (15)) by explicitly filtering the resolved solution on the 643 grid. We use a Gaussian filter with a width of 4∆f , where ∆f refers to the grid spacing on the 643 grid. (2) The similarity model is evaluated by explicitly filtering the fields on the 323 grid. We use a Gaussian filter with a width of 2∆c , where ∆c refers to the grid spacing on the 323 grid. (3) The nonlinear gradient model is evaluated using the fields on the 323 grid.
4.1. Estimates for the SGS terms In this subsection, we evaluate the SGS terms using a snapshot from the direct numerical simulation at t = 1.9. We compare the predicted SGS terms with an estimate obtained from the fully resolved calculation, and with
Although the difference in resolution between the coarse and fine grid is small, it is sufficient for permit an initial test of the SGS models. Successful models at the current resolution do not guarantee good performance in subsequent calculations when the resolution of the
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
115
Fig. 8. Intensity of divergence of the SGS heat flux (upper panels) and z-component of divergence of the SGS momentum flux on the mid-depth of the simulation domain (e.g. z = 0). The result obtained by direct estimation on the resolved grid is given in the left, that obtained by the similarity model is plotted in the middle, and the result of the gradient model is plotted in the right panel.
coarse and fine grids are substantially different. On the other hand, failure at the current resolution would almost certainly ensure failure at more disparate resolutions. The results for the divergence of the SGS heat flux −∇ · I and the z-component of the other SGS terms −∇T , E−1 Pm−1 ∇M, and a at the mid-depth of the domain are shown in Figs. 8 and 9. As seen in Figs. 8 and 9, general features of these three predictions are similar to each other. These comparisons indicate that the nonlinear gradient model is nearly equivalent to the similarity model, and that both models provide a good description of the direct estimation in the mid-plane of the simulation domain. To compare these predictions in more detail, we plot the spectra of the SGS terms with respect to the ydirection in Fig. 10. We can find many common characteristics among all of the SGS terms obtained by the nonlinear gradient model. For example, the SGS terms of the gradient model have larger small-scale components than the predictions of the similarity model. Because of these differences, the small-scale components of the nonlinear gradient model are in better agreement with the small-scale components of the direct estimation.
Only the SGS Lorentz term deviates substantially from the direct estimation. The SGS Lorentz term predicted by the similarity model is even worse. These discrepancies could be reduced by changing the magnitude Csim for the SGS Lorentz term. In general, the magnitude of each SGS term is approximately 10% of the resolved nonlinear terms in the governing equations. One exception occurs in the SGS heat flux, where the SGS term is approximately 4% of the resolved heat flux. The reason of this difference is that the spectrum of the resolved heat flux has a steep slope compared with the spectra of the other resolved terms. We also note that SGS Lorentz force has much larger amplitude than the Reynolds stress. The SGS terms predicted around boundaries are different from those predicted in the interior. The intensity of the divergence of the SGS heat flux and z-component of the SGS Lorentz term on a surface near the lower boundary (z = −0.4975) at t = 1.9 is shown in Fig. 11. The pattern of z-component of SGS Lorentz term is similar to the direct estimation, but the intensity of the SGS terms are smaller than the direct estimates. The divergence of heat flux from the nonlinear gradient model
116
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
Fig. 9. Intensity of z-component of SGS Lorentz term (upper panels) and z-component of SGS induction term (lower panels) term on the mid-depth of the simulation domain (e.g. z = 0). The result obtained by direct estimation on the resolved grid is given in the left, that obtained by the similarity model is plotted in the middle, and the result of the gradient model is plotted in the right panel.
differs in pattern and amplitude from the results of both the direct estimation and similarity model. These differences are most evident in the top-right corner of the plots, where intense negative values are predicted by the nonlinear gradient model, whereas intense positive values are obtained in the direct estimation and similarity model. These results suggest that the SGS terms in the nonlinear gradient model require modification around boundaries. 4.2. Large-eddy simulation using the SGS model As described in the previous section, the SGS terms can be reliably approximated using the nonlinear gradient model in the region away from the boundary layer. To test the utility of this model, we perform a large-eddy simulations with a 323 grid from t = 1.9 to 2.5, using the snapshot at t = 1.9 as an initial value. We compare results obtained from three simulations: (1) A direct simulation with 643 nodes (resolved DNS). (2) A direct simulation with 323 nodes (unresolved DNS).
(3) A simulation using the nonlinear gradient model with 323 nodes (LES). We plot time evolution of the kinetic and magnetic energy averaged over simulation domain for the three cases in Fig. 12. A comparison of the magnetic energy in Fig. 12 shows that the resolved DNS and LES are quite similar prior to t = 2.05. Both of these solutions have large increases in the magnetic energy at t = 2.05, whereas the increase in the unresolved DNS is much smaller. However, the behavior of the LES deviates from the resolved DNS case after t = 2.1. Because of the large time variation in the kinetic and magnetic energies, we compare the time average of the energies from t = 2.2 to 2.5 to avoid the effects of the initial conditions (see Table 1). The average kinetic and magnetic energies in the LES is intermediate between the resolved DNS case and the unresolved DNS case. To compare the convection and magnetic field pattern, we choose a cross section at the mid-plane of the domain after taking 20,000 time steps from the initial value. The corresponding solution time is t = 1.95, which is prior to the divergence of the magnetic energy at t ∼ 2.05. As
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
117
Fig. 10. One-dimensional spectra of SGS terms in the mid-depth of the domain (z = 0.0) at t = 1.9. The prediction of the nonlinear gradient model is plotted by circles with the thick line, while the result of the similarity model is plotted by triangles. The direct estimate on the 643 grid is plotted by ‘+’, and the corresponding resolved terms are plotted by cross with dotted line for reference.
Fig. 11. Intensity of divergence of the SGS heat flux (upper panels) and z-component of the SGS Lorentz term (lower panels) on a surface near the lower boundary (e.g. z = −0.4975). The direct estimation is given in the left, while the result of the similarity model is plotted in the middle. The result of the gradient model is plotted in the right panels.
118
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123 Table 1 Average kinetic and magnetic energies for 2.2 < t < 2.5
Fig. 12. Time evolution of the kinetic energy (top panel) and the magnetic energy (bottom panel) for the resolved DNS, unresolved DNS, and the LES using the nonlinear gradient model.
Run
Resolved DNS
Unresolved DNS
With SGS model
Ekin Emag
7500.3 42922.2
8070.1 49139.3
7727.0 43382.4
seen in Fig. 13, we can find similarities in the pattern of the fields, although the details differ. In general, we can see that the unresolved DNS case has more power in small scale fields. To quantitatively investigate the similarity of the field patterns, we plot the time evolution of the cross correla¯ z in the mid-plane with the resolved DNS tion of the B results in Fig. 14. The correlation coefficient decreases rapidly to 0 by t = 1.95 in both cases, and it is difficult to distinguish between the LES result and unresolved DNS result on the basis of this diagnostic (see Fig. 14). On the other hand, the phase angle between the average horizontal magnetic field and the x-axis is almost same for each model at t = 1.95 (see Fig. 15). These results suggests that the small scale components deviate rapidly because the correlation is strongly influenced by the small scale fields. However, the phase angle also deviates after t = 2.05 (see Fig. 15). The phase angle of the
Fig. 13. Intensity of the z-component of vorticity (upper panels), and magnetic field (lower panels) in the mid-plane of the domain (z = 0.0) at t = 1.95. The result obtained by the resolved DNS is given in the left, the result obtained by the unresolved DNS is plotted in the middle, and the result by the LES is plotted in the right panel.
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
119
Fig. 14. Time evolution of the correlation coefficient for the zcomponent of the magnetic field in the mid-plane. The correlation between the resolved and unresolved DNS is plotted by dotted line, and the correlation between the resolved DNS and LES is plotted by solid line.
LES rotates more slowly than that in the resolved DNS case, whereas the phase angle in unresolved DNS case changes faster than that in the resolved DNS result. To investigate the effects of the SGS terms on the resolved fields, we plot a power spectrum of the kinetic and magnetic energies with respect to the y-direction for t = 2.4 in Fig. 16. The spectrum of the kinetic energy in the LES is much closer to the kinetic energy spectrum in resolved direct simulation. The unresolved DNS substantially overestimates the kinetic energy at the small scales. This result suggests that the SGS terms provide adequate dissipation of the kinetic energy at the small scales. The spectra of magnetic energy have similar characteristics, but the small scale components of the magnetic energy in the LES become much smaller than those in the resolved DNS case. These characteristics also can be seen the snapshot of the fields at t = 1.95 (see Fig. 13). Recalling the average kinetic and magnetic energies in Table 1, the difference of the averaged energy is due mainly to the difference in the large scale fields. The relative importance of the SGS terms in the momentum equation can be inferred from the power spectra
Fig. 15. Time evolution of the phase angle of the average horizontal magnetic field for the resolved DNS (thick solid line), unresolved DNS (dashed line), and the LES using the nonlinear gradient model (thin solid line).
Fig. 16. One-dimensional spectra of the kinetic and magnetic energies at t = 2.4.
of the SGS terms (see Fig. 10). The SGS Lorentz term is approximately 10 times larger than the Reynolds stress. In fact, the SGS Lorentz term is comparable to the resolved inertial term. The magnitudes of the SGS Lorentz and Reynolds stress suggest that the Reynolds stress has a smaller role in the present simulation, and that the difference in the kinetic energy spectrum between the LES and unresolved DNS cases is due to the SGS Lorentz term. Energy balances are also important for understanding the role of the SGS terms. The conservation of the resolved kinetic energy is obtained by integrating the product of u¯ and Eq. (3) over the fluid volume. The resolved magnetic energy can similarly be described by ¯ and Eq. (5) over the fluid integrating the product of B volume. The result is ∂ 1 2 u¯ dV ∂t 2 1 2 =− u¯ u¯ + u¯ χ¯ + ω¯ × u¯ · dS 2 ¯ + −ω¯ 2 + RaE−1 u¯ z (T¯ − T0 ) + Pm−1 E−1 u· ¯ − u¯ · ∇Tij + Pm−1 E−1 u¯ · ∇Mij dV, (J¯ × B) (40)
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
120
Fig. 17. Time evolution of energy flux for the kinetic (upper panel) and magnetic energies (lower panel) for the LES. It is noted that all of the SGS contributions are negative. The work done by the Reynolds stress is not shown in the upper panel because it is not distinguishable from zero in this plot.
∂ ∂t
1 ¯2 B dV 2
=−
¯ ×B ¯ dS E
[−Pm−1 J¯ 2
¯ − J¯ · a] dV, − J¯ · (u¯ × B)
(41)
Fig. 18. Schematic illustration of the role of the SGS induction and Lorentz terms in the kinetic and magnetic energies. The work done by the resolved Lorentz force transfers kinetic energy from the resolved velocity into the resolved magnetic energy. Similarly, the SGS Lorentz term draws kinetic energy from the resolved velocity into magnetic energy of the unresolved field. Generation of magnetic energy by the resolved fields is positive, whereas the SGS induction term dissipates magnetic energy.
dicating a sink of resolved magnetic energy. By contrast the resolved induction term is positive, indicating a constructive contribution to the magnetic energy. We summarize the energy transfer between resolved and unresolved component by the SGS terms in Fig. 18. The SGS Lorentz term and SGS induction term transfer energy from the resolved components of the kinetic and magnetic energies into the unresolved components.
where ¯ = Pm−1 J¯ − (u¯ × B) ¯ − a. E
(42)
The boundary conditions (Eqs. (30)–(33)) ensure that the surface integrals in Eqs. (40) and (41) vanish in the present model. The importance of the various terms in Eqs. (40) and (41) is shown in Fig. 17. The SGS Lorentz force and induction term have significant role in the energy balance because they are large compared with the other SGS terms. As seen in Fig. 17, the work done by the resolved and SGS Lorentz force has same sign through the simulation. On the other hand, the SGS induction term is negative throughout the simulation in-
5. Discussion The results described in the previous sections indicate that the SGS model transfers energy from large to small scales. The SGS Lorentz term is primarily responsible for dissipating kinetic energy in the resolved scales, whereas the SGS induction term dissipates large-scale magnetic energy. Comparison of the magnetic energy spectra in Fig. 16 suggests that the SGS induction term draws too much energy from the resolved components. The fact that the nonlinear gradient model reproduces the spectral amplitude of the direct estimate in Fig. 10 suggests
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
121
Fig. 19. Intensity of the y-component (upper panels) and z-component (lower panels) of the SGS heat flux on a surface near the lower boundary (e.g. z = −0.4975). The direct estimation is given in the left, while the result of the similarity model is plotted in the middle. The result of the gradient model is plotted in the right panels.
that the phase of those two terms is slightly different. One way to reduce the dissipation caused by the SGS induction term is to reduce the model constant Csim in Eq. (27), although the spectral amplitudes of the SGS induction term would then be small. Similar problems arise with the SGS Lorentz term. The amplitude of the SGS Lorentz term is much smaller than that predicted by the direct estimation (see Fig. 10), yet the dissipation of the kinetic energy is about right. This inconsistency may be a consequence of comparing spectra at a specific snapshot in time. The spectra of the SGS terms relate to the time derivative of the energies, rather than the energies directly. It would be preferable to base comparisons on time-averaged spectra. Another problem with the current implementation of the nonlinear gradient model occurs near the boundaries. We previously noted in Fig. 11 that the predicted pattern of −∇ · I near the boundary differs from that of the direct estimation. In order to identify the origin of this discrepancy, we plot the components Iy and Iz near the boundary at the same time step. We can see that these components are in much better agreement with the SGS heat flux from the direct estimation, although the
amplitudes are slightly different. The nonlinear gradient model appears to overestimate Iy and underestimate Iz . If we repeat this comparison at different distances from the boundaries, we find that the discrepancy in the amplitude of the heat flux components varies systematically with distance. Specifically, the nonlinear gradient model increasingly overestimates the direct estimate as we approach the boundary. This suggests that the model coefficients Csim should decrease toward the boundary. Germano et al. (1991) devised a scheme to obtain the model coefficient Csim by explicitly filtering the resolved fields. This scheme was originally developed to estimate the coefficient for the Smagorinsky model. The success of the so-called dynamic model was demonstrated by Piomelli (1993), who showed that the coefficient of the Smagorinsky model approaches to 0 toward the boundary. The dynamic model is quite general and can be applied to other SGS models, including the nonlinear gradient model. Thus, the potential exists to evaluate the coefficients of the nonlinear gradient model using a dynamic scheme. If this approach can be used to approximate the constant Csim (z) that we obtained from a fit to the direct estimation of the SGS terms, then we
122
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
Fig. 20. Comparison of the SGS heat flux when calculated using −∇ · (uT ) + ∇ · (u¯ T¯ ) (left panel) and −[∇ · (uT − u¯ T¯ )] (middle panel). Both terms are evaluated on the 643 grid near the lower boundary (i.e. z = −0.4975). The result of the nonlinear gradient model on the 323 grid is plotted in the right panel for comparison.
can expect good performance in large-eddy simulations of convection-driven dynamos. Another source of difficulty near the boundaries may be caused by the use of a non-uniform grid. We have seen in Fig. 19 that the pattern of Iy and Iz of the gradient model around the boundary is quite similar to that of the direct estimation and the similarity model. However, the pattern of −∇ · I by the gradient model is quite different from that obtained by the direct estimation and similarity model (see Fig. 11). Ghosal and Moin (1995) pointed out that the filtering operation does not, in general, commute with the operation of the differentiation if we use a filter with a variable filter width. Commutation of the operations is customarily assumed in deriving the governing equations for the large-scale fields. For example, in the heat equation for the large scale temperature field, we assume that ∂ ∂ (ui T ) = (ui T ). ∂xi ∂xi
(43)
In the present study, the size of element changes with respect to the horizontal position (see Fig. 2) and filter width is changed by the size of the element in the physical space. To asses the error caused by interchanging the order of differentiation and filtering, we compare two representations of the SGS heat flow. One representation is given by −∇ · (uT ) + ∇ · (u¯ T¯ ), which was previously used to define the direct estimate in Fig. 11. The other representation is given by −[∇ · (uT − u¯ T¯ )]. Both terms are evaluate near the boundary, where the grid spacing is most uneven. Differences between the two terms provide a measure of the commutation error Fig. 20. We note that the pattern of −[∇ · (uT − u¯ T¯ )] is similar to −∇ · I sim , as obtained using the nonlinear gradient model. Consequently the error in the SGS model near the boundaries
may be a result of the commutation error. Vasilyev et al. proposed an approach to construct a commutative filter for the unstructured grid. Construction of a commutative filter will be important for improving the SGS model especially near the boundaries. 6. Conclusion We implemented a sub-grid scale model (SGS) in a dynamo simulation code using the finite-element method. We applied the nonlinear gradient model which is a version of the similarity SGS model to the four nonlinear terms in the MHD equations, e.g. the heat and momentum fluxes, Lorentz term, and the induction terms. First, we test the predictions of the SGS model using snapshots from a direct numerical simulation (DNS). The result showed that the spatial pattern of the SGS terms in the nonlinear gradient model is similar to that obtained by the direct estimation from a resolved DNS at the mid-plane. The only exception occurs with SGS Lorentz term, where the nonlinear gradient model underestimates the direct estimation. Simulations using this SGS model are performed in a rotating plane layer model, and we compared the results with the output a direct numerical simulation using either the same or finer spatial resolution. All of those simulations have a large time variation in the kinetic and magnetic energies. However, the time averaged kinetic and magnetic energies energy obtained by LES are intermediate between the resolved DNS case and the unresolved DNS case. Comparing field patterns, the patterns of convection and magnetic field in the LES are more like those in the resolved DNS case than those in the unresolved DNS because the unresolved results have strong small-scale fields.
H. Matsui, B.A. Buffett / Physics of the Earth and Planetary Interiors 153 (2005) 108–123
We also investigate the energetics of the convection and magnetic field in the simulation. The work done by the SGS Lorentz term and the magnetic energy generated by the SGS induced terms are both negative, indicating that energy is transferred form the resolved components to the unresolved components. While the simulations suggests that the SGS models works adequately, we also find some discrepancy between the LES and resolved DNS results. The spectrum of the kinetic energy in the LES is much closer to that in the resolved DNS (as compared with the unresolved DNS). On the other hand, small-scale components in the magnetic energy of the LES are smaller than those in the both resolved and unresolved DNS, implying that the SGS model is overly dissipative. The largest discrepancies in the SGS models occur near the boundaries. To improve the present model, we need to modify the filter to reduce the commutation error caused by interchanging the order of filtering and spatial differentiation. The dynamic model also has the potential to improve the prediction of SGS terms around the boundaries by applying the model coefficients as a function of the vertical position z. Acknowledgements The work is supported by the National Science Foundation. The authors would like to the National institute for Polar Research and to Center for Planning and Information systems in Japan Aerospace Exploration Agency for performing the present simulations. References Braginsky, S.I., Meytlis, V.P., 1990. Local turbulence in the earth’s core. Geophys. Astrophys. Fluid Dyn. 55, 71–87. 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. Buffett, B.A., 2003. A comparison of subgrid-scale models for largeeddy simulations of convection in the earth’s core. Geophys. J. Int. 153, 753–765. Christensen, U.R., Olson, P., Glatzmaier, G.A., 1999. Numerical modelling of the geodynamo: a systematic parameter study. Geophys. J. Int. 138, 393–409. Germano, M., 1986. A proposal for a redefinition of the turbulent stresses in the filtered navier stokes equations. Phys. Fluids 29, 2323–2324. Germano, M., Piomelli, U., Moin, P., Cabot, W., 1991. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids 3, 1760–1765. Ghosal, S., Moin, P., 1995. The basic equations for the large eddy simulation of turbulent flows in complex geometry. J. Comput. Phys. 118, 24–37.
123
Glatzmaier, G.A., Roberts, P.H., 1995. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle. Phys. Earth Planet. Interiors 91, 63–75. Glatzmaier, G.A., Roberts, P.H., 1997. Simulating the geodynamo. Contemp. Phys. 38, 269–288. Jones, C.A., Roberts, P.H., 2000. Convection-driven dynamos in a rotating plane layer. J. Fluid Mech. 404, 311–343. Kageyama, A., Sato, T., the Complexity Simulation Group, 1995. Computer simulation of a magnetohydrodynamic dynamo. Phys. Plasmas 2, 1421–1431. Kuang, W., Bloxham, J., 1997. An earth-like numerical dynamo model. Nature 389, 371–374. Kuang, W., Bloxham, J., 1999. Numerical modeling of magnetohydrodynamic convection in a rapidly rotating spherical shell: weak and strong field dynamo action. J. Comput. Phys. 153, 51–81. Leonard, A., 1974. Energy cascade in large-eddy simulations of turbulent fluid flows. Adv. Geophys. 18, 237–248. Lesieur, M., Metais, O., 1996. New trends in large-eddy simulations of turbulence. Annu. Rev. Fluid Mech. 28, 45–82. Matsui, H., Okuda, H., 2004. Development of a simulation code for mhd dynamo processes using the geofem platform. Int. J. Comput. Fluid Dyn. 18 (4), 323–332. Matsushima, M., 2001. Expression of turbulent heat flux in the earth’s core in terms of a second moment closure model. Phys. Earth Planet. Interiors 104, 137–148. Matsushima, M., Nakajima, T., Roberts, P.H., 1999. The anisotoropy of local turbulence in the earth’s core. Earth Planets Space 51, 277–286. Meneveau, C., Katz, J., 2000. Scale-invariance and turbulence models for large-eddy simulation. Annu. Rev. Fluid Mech. 32, 1–32. Moffatt, H.K., Proctor, M.R.E., 1982. The role of the helicity spectrum function in turbulent dynamo theory. Geophys. Astrophys. Fluid Dyn. 21, 265–283. Parker, E.N., 1955. Hydromagnetic dynamo models. Astrophys. J. 122, 293–314. Phillips, C.G., Ivers, D.J., 2001. Special interactions of rapidly-rotating anisotropic turbulent viscous and thermal diffusion in the earth’s core. Phys. Earth Planet. Interiors 128, 93–107. Phillips, C.G., Ivers, D.J., 2003. Strong field anisotropic diffusion models for the earth’s core. Phys. Earth Planet. Interiors 140, 13–28. Piomelli, U., 1993. High reynolds number calculations using the dynamic subgrid-scale stress model. Phys. Fluids A 5, 1484– 1490. R¨adler, K.-H., Wiedemann, E., Brandenburg, A., Meinel, R., Tuominen, I., 1990. Nonlinear mean-field dynamo models: stability and evolution of three-dimensional magnetic field configurations. Astron. Astrophys. 239, 413–423. Sakuraba, A., Kono, M., 1999. Effects of the inner core on the numerical simulation of the magnetohydrodynamic dynamo. Phys. Earth Planet. Interiors 111, 105–121. Shimizu, H., Loper, D.E., 1997. Time and length scales of buoyancydriven flow structures in a rotating hydromagnetic fluid. Phys. Earth Planet. Interiors 104, 307–329. Winckelmans, G.S., Wray, A.S., Vasilyev, O.V., Jeanmart, H., 2001. Explicit-filtering large-eddy simulation using the tensor-diffusivity model supplemented by a dynamic smagorinsky term. Phys. Fluids 13, 1385–1403. Zhang, K.K., Jones, C.A., 1997. The effect of hyperviscosity on geodynamo models. Geophys. Res. Lett. 24 (22), 2869–2872.