2. Model and Method of Numerical Calculation |

A 2-dimensional non-divergent flow on a rotating sphere governed by the following vorticity equation will be examined.

Here,

ψ(λ, μ, t): stream functions;ω(λ, μ, t): vertical component of vorticity (≡ Δψ); λ: longitude; μ: sine latitude; t: time;Ω: angular speed of rotation of sphere;ν: hyper-viscosity coefficient, and_{2p}p: order of hyper-viscosity. Furthermore, Δ is the horizontal Laplacian and is defined as follows.Note that,

- the dimensionless units of length and time are the radius of the sphere and a given time scale T, respectively,
- the above vorticity equation corresponds to the Lagrangian conservation law for potential vorticity of
q≡ω+ 2Ωμ in cases of inviscid fluids,- the hyper-viscosity coefficient term on the right-hand side of the vorticity equation has been adopted for convenience to prevent numerical instability and failure, and
- the "+2" in the hyper-viscosity coefficient term on the right-hand side of the vorticity equation (Δ+2) has been added to conserve total angular momentum.
The dependent variable

ψis expanded using spherical harmonics aswhere,

P(μ) is the associated Legendre functions normalized by 2, and_{n}^{m}- the reason the discretization starts from n = 2 is that the component of n = 1 corresponds to the total angular momentum, which is the conserved quantity in the vorticity equation.
Based on the above expansion, the simultaneous ordinary differential equations for the expansion coefficient

ψ(t) is yielded from the vorticity equation using the spectral Galerkin method. The non-linear terms are calculated in real space, and a technique is adopted for returning the results back into the wavenumber space by Fourier-Legendre transformation. To eliminate aliasing, global surface domain is covered by a transform grid of 2,048 (longitudinal) × 1,024 (latitudinal) when the truncation wavenumber N = 682, and a grid of 1,024×512 when N = 341._{n}^{m}In the experiment for decaying turbulence, an initial turbulent field will be given to observe time evolution. The energy density is given by

The initial energy spectrum is constrained to be the distribution function employed by Cho and Polvani (1996)

and each band wavenumber component

mwill be given random amplitude and phase. Here, coefficientAwill be set so that the total dimensionless kinetic energyE≡ Σ_{total}_{n=2}^{N}E(n,t=0) will be equal to 1. n_{0}is a parameter that gives the central wavenumber, and will take the values of either 10, 50, or 100. γ is a parameter that determines the spectral width, and a larger value for γ will confine the initial energy to a narrower band. In the present study, γ was set to 1,000.

- The fact that
E= 1 is adopted conversely implies that (radius of sphere) / √_{total}Ehas been adopted as the time scale for making the vorticity equation dimensionless._{total}The angular velocity of rotation

Ωwill be used as the experimental parameter and given six values, 0, 25, 50, 100, 200, and 400. WhenΩ≠ 0, the planet will makeΩ/(2π) rotations during unit time.

- In this scaling scheme, the relationship between dimensionless
Ωand dimensional spherical radius ofa, angular velocity of rotation of^{*}Ω^{*}, and mean wind velocity ofUcan be represented as^{*}Ω= √ 2 ×a^{*}Ω^{*}/U. Since the model used in the present study is a 2-dimensional incompressible fluid, it is far from being a precise representation of the actual planetary atmosphere, but if forced to make analogies with actual conditions, the above values for Earth are^{*}a= 6.37×10^{*}^{6}m,Ω^{*}= 7.29×10^{-5}s^{-1}, and if 15 ms^{-1}is chosen for the value ofU, then^{*}Ω= 43.8. For Jupiter, these values area= 6.99×10^{*}^{7}m,Ω^{*}= 1.76×10^{-4}s^{-1}, and if 50 ms^{-1}is chosen for the value ofU, then^{*}Ω= 348. Therefore, if we were to choose the results corresponding to the cases for Earth and Jupiter, then the experiments forΩ= 50 andΩ= 400 will roughly have parameters corresponding to that of Earth and Jupiter, respectively.Note here that the wavenumber at which the non-linear term of the vorticity equation and the linear rotational effect term are nearly equal (wavenumber corresponding to the Rhines scale), n

_{β}≡ √ (πΩ/(4U))= √ (πΩ/(4√(2×E)) ) =√(π_{total}Ω/(4√ 2)), will be an important indicator. Actually, the manner of time evolution of the decaying turbulence is strongly dependent on the magnitude relation between n_{β}and the initial energy peak wavenumber n_{0}. When n_{β}> n_{0}, then the linear term will dominate over the non-linear term from the initial stages, and so the inverse energy cascade predicted by the turbulence theory will not always appear.For the numerical viscosity term,

pis set to 8, andν= 1×10_{2p}^{-43}(for N = 682), and 1×10^{-38}(for N = 341). Time integration is done using the 4th-order Runge-Kutta method, with the time increment being set to Δt = 5×10^{-4}(for N = 682) and 1×10^{-3}(for N = 341), and time evolution will be computed to t = 5.For the numerical calculations, the efficient transform code of ISPACK (ishioka, 1999) has been employed.

<<index <previous next>