Multiple space and time analysis of the barotropic equation

Florin Spineanu Cadarache 2001 and Magurele 2002

Abstract:This part is based on the derivation provided by Tan and Tan&Boyd (B. Tan, J. Atmospheric Sciences 53, 1604 (1996), B. Tan and J. P. Boyd, Chaos, Solitons & Fractals 11, 1112 (2000) ). We extend the calculation by joining numerical realization.

Multiple time and space scale analysis of the ion mode equation reduced to the barotropic form

The quasigeostrophic barotropic atmospheric model leads to the barotropic equation (see Horton, Phys.Reports where the quasigeostrophic approximation is explained in relation with the Ertel's theorem). This equation has been studied in the atmosphere science by means of the multiple space and time scales method. Our aim is similar, i.e. to obtain an equation for the envelope of the fluctuating field $\phi $ in order to study the possible poloidally nonuniform profiles of the turbulence averaged amplitude. The connection between the original equation and the final slow-time and large-spatial scales equations is not simple : much numerical work is required in order to connect the input physical parameters of the plasma with the formal coefficients in the final equations. The detailed analytical work is presented in Ref(Tan) and the necessary formulas are reproduced here for convenience. To simplify the expressions the following notation is introduced MATH and the equation can be written MATH

The analysis on multiple scales starts by introducing the variables on scales separated by the parameter $\varepsilon $: MATH MATH This gives for the derivatives MATH MATH

The solution is adopted to be of the form MATH We denote the linear part of the operator in the equation by MATH Substituting the Eq.(psiserie) and the forms of the derivation operators Eqs.(dertserie, derxserie) in the equation of motion we obtain from the equality between the coefficients of the powers of the variable $\varepsilon $: MATH MATH MATH

The solution is represented on the space and time scales.

The enveloppe equation

The solution of the linear part of the operator, Eq.(linear) is adopted as the superposition of two propagating waves MATH The two functions $\varphi _{n}$, $n=1,2$, are solutions of the equation MATH with MATH In this formula, $\left( 0,L\right) $ are the limits of the spatial domain in the minor radius direction where the flow MATH ; MATH are the phase velocities. To find the solutions of this equation we have to use numerical methods. This is a Schrodinger-type equation where the potential depends on the energy, (since the phase velocity $c_{n}$ depends on the wavenumber $k_{n}$). We must use numerical methods suitable for Sturm-Liouville equations and in general the solution exists for only particular values of the parameters. We choose to fix the frequencies $\omega _{n}$ and consider $k_{n}$ as the parameter to be determined. The problem is complicated by the possible occurence of singularities in the potential due to resonances between the phase velocity and the fixed zonal flow. We will avoid any resonance since the physical content of the processes consisting of the direct energy transfer between the fluctuations and the flow represents a different physical problem and requires separate consideration. In the present work we examine the effect of the turbulence on the poloidal rotation as being due to the nonuniform diffusion rates and Stringer mechanism.

The parameters in the Eqs.(phin) are $\beta $, $\varepsilon $ and the sheared flow velocity $U(x)$. We take for $U\left( x\right) $ a symmetric Gaussian-like profile shifted in amplitude and retain its maximum, half-width and asymptotic value (i.e. at $x=0$ or $L$) as parameters. The two frequencies have fixed values during the eigenvalue search but they must be considered free parameters as well. This builds up a large space of parameters which should be sampled. For some point in this space the eigenvalue problems Eqs.(phin) are solved and the phase velocities compared with the profile of $U\left( x\right) $. The presence of resonances renders the set of parameters useless. Since (as will result from the formuals below) there are also other possible resonances, the determination of a useful solution MATH is difficult and requires many trials.

The numerical methods that have been used in solving Eqs.(phin) belong to two classes: boundary value integration methods and respectively shooting methods. When the sign of the ``potential'' in the equation is negative everywhere on $(0,L)$ , which occurs for short wavelengths and smooth fluid velocity profile the boundary conditions cannot be fulfilled except for asymtotically, i.e. assuming arbitrary small nonzero values at these limits, since the solution decays exponentially. However most of the boundary value integrators we have tried gave amplitudes of magnitudes comparable to the boundary value, which essentially are homotopic to the trivial vanishing solution (the equation is homogeneous) and render useless this approach. We have to use a shooting method accepting the difficulty that the initial value of the parameter to be determined (eigenvalue) must be placed not far from the final value. We have used the specialized package SLEIGN2 interactively. For the investigation of ranges of parameters most of the calculations have been performed using the routine D02KEF of NAG. For identical conditions the two numerical codes gave the same results within the accepted tolerance.

Using the routine D02KEF we divide the integration range $\left( 0,L\right) $ into four intervals corresponding roughly to a different behaviour of the potential. This should simplify the task of integrator. A single matching point is assumed, usually at the centre. The tolerence is $10^{-8}$ which however does not exclude fake convergence in case of an initialisation of the eigenvalue very far from the true solution..

We notice from Eqs.(phin) that the potential can easily be dominated by high shear $U^{\prime \prime }$ especially for smaller wavenumbers $k_{n}$. Actually we are more interested to investigate regimes where only mild effects of the seed sheared rotation arise during the process of self-modulation since, as explained before, the possible evolution to the nonuniformity in the poloidal direction will be the source of the torque on the plasma. Avoiding regimes of $U^{\prime \prime }$-dominated potentials we assume low shear and take the Gaussian profile with MATH ($m/s$) represented in Figure (figu).
u_18.ps
Profile of the seed sheared poloidal velocity. Also are plotted the first and the second derivatives.

The following set of physical parameters has been used: minor radius $a=1\;m$, magnetic field $B_{T}=3.5\;T$, electron temperature $T_{e}=1\;KeV$, ion temperature $T_{i}=1\;KeV$, density $n=10^{20}\;m^{-3}$, $L_{n}=10\;m$. The parameters in the barotropic equation are then: $\ \beta =33.526$ and MATH. The two frequencies are MATH and MATH. It is assumed that the space region of significant magnitude of $U\left( x\right) $ is $L=0.1\;m$. The following two eigenvalues are obtained, normalised to $L^{-1}$: $k_{1}=0.5525$ and $k_{2}=3.2029$. The potentials in the Schrodinger-like equations and the eigenfunctions are shown in Figure (figphi).
phi_18.ps
The potentials in the Schrodinger-type equations and the eigenfunctions $\varphi _{1,2} $. The derivatives of the eigenfunctions are aslo plotted.

The two eigenfunctions $\varphi _{1,2}$ have an important role in the following steps of the calculation. Substituting the Eq.(psiat1) in the barotropic equation and expanding the operators, we obtain: MATH

where MATH We have $k_{12}=3.755$, $\omega =8.585$. MATH or: MATH, MATH. MATH MATH MATH MATH These functions can be obtained after calculating the eigenfunctions $\varphi _{1,2}$ and are represented in Figure (figg).
g_18.ps
Graphs of the functions $g$'s.

Various solutions to the equation (lpsi2) are possible, taking the space-time dependence similar to one of the terms which compose the right hand side.

The second term suggests solutions of the type MATH

The third term suggets solutions of the type MATH

The fourth term gives MATH

In the formulas above the functions MATH are solutions of the equations: MATH MATH MATH These equations are integrated numerically with a boundary value integrator, a finite difference method. The results are shown in Figures (figw1), (figw2) and (figw3). In the following figures are explained also the functions $W_{1n}$ which will be introduced in Eq.(eqw1n).
y1_18.ps
Graphs of the potentials occuring in the equations for the functions $w$'s.


y2_18.ps
Graphs of the source terms in the equations for the functions $w$'s.


y3_18.ps
Graphs of the functions $w$'s.

The first term gives a solution of the following form MATH where the function MATH satisfies the equation MATH

This equation serves to generate a condition on the two amplitudes $A_{1,2}$.

When the left hand side of the equation is multiplied by the function $\varphi _{n}$ and integrated between the limits on the $x$ domain, it gives zero, due to the boundary conditions. The same must then be true for the right hand side. This gives a ``solubility condition'', on the space time scales slower at order $1$: MATH where MATH Expressed in units of length $L$ and units of time $t_{0}=10^{-4}$, the group velocities are : $c_{g1}=0.30328$ and $c_{g2}=3.9908$.

Then it is clear that the amplitude $A_{n}$ propagates at the speed $c_{gn}$ : MATH Now we introduce a new function MATH The equation satisfied by $W_{1n}$ is MATH with the bounday conditions MATH

At this order $\left( 2\right) $ there are the following solutions MATH The last function will be determined by the condition of solubility at the next order on the space-time scales (MATH). To obtain the evolution equations for $\Phi $ and $A_{n}$, we use the expressions for MATH Eq.(psiat1) and MATH, Eq.(psiat2) in the right hand side of the equation (lpsi3), and all inhomogeneous terms are obtained. In additions there are terms which are not dependent on $y$ and $t$, and these terms should be zero. This leads to the solubility condition relating the function $\Phi $ (correction to the mean flow due to the presence of the wave packets) to the amplitude $A_{n}$. MATH This relation will be used later.

Returning to the equation (lpsi3), we examine the content of the inhomogeneous part (the right hand side). There are terms which have a space-time dependence MATH and in order to have solutions of this type, we need to impose solubility conditions, otherwise there will arise secular terms at the order MATH. To find these restrictions (solubility conditions) we multiply the equation with the function conjugated to (expkn) MATH and integrate time on the interval MATH and space: $y$ on the interval MATH and $x$ on the interval $\left( 0,L\right) $. On the left hand side we obtain zero, so this is what we must obtain also from the right hand side. The conditions which results are: MATH MATH The notations which have been introduced are: MATH, MATH, MATH, MATH, MATH and MATH MATH MATH

The equations satisfied by all these functions are given given in Ref.([Tan]). MATH MATH MATH MATH MATH MATH MATH We use a boundary value integrator for the inhomogeneous equations and (the NAG routine D02JAF). The results are given in the Figures (figf1) and (figf2).
f1_18.ps
Graphs of the first group of functions $f$.


f2_18.ps
Graphs of the second group of functions $f$.

After numerical integrations the follwing values are obtained for the constants: MATH, MATH, $\rho _{1}=-0.10182$, $\rho _{2}=1.8115$, MATH, MATH.

The following transformation (Jeffrey and Kawahara) is performed in Ref.( [Tan]) : MATH Then the equation of connection between $\Phi $ and $A_{n}$ becomes MATH

The solution of the equation (eqphia) is considered to be of the form MATH and the equation for MATH results: MATH

The change of variables is performed in the equations for the amplitudes $A_{n}$ and, again, the Eq.(phimarea) is used: MATH MATH where MATH MATH MATH MATH

The following transformation is convenient: MATH This leads to the following equations: MATH MATH

The constants are: MATH, $\sigma _{2}=1.9331$, $\nu _{12}=-0.53387$, $\nu _{21}=-2.9014$. The constants $\alpha _{1,2}$ are called the ``dispersion'', $\sigma _{1,2}$ are called the nonlinearity (``Landau'') constants ; $\nu _{12,21}$ represent the coupling of the two amplitudes.

Numerical experience shows that the dispersion constants are sensitive to the precision of the integration scheme. We use for D02KEF the value for tolerence $10^{-8}$.

The system of coupled Cubic Nonlinear Schrodinger Equations has been integrated numerically in relation with many applications. One important conclusion is that, in the case of a single equation, the sign of the product of the dispersion and of the nonlinearity parameters establishes the possibility of a soliton formation.

Numerical search for classes of admissible envelope equations

In principle an interactive search for eigenvalues (either with SLEIGN or D02KEF) is able to obtain the eigenvalues $k_{1},k_{2}$ eventually after series of less-inspired initializations. However in many situations resonances occur and the sources and/or the potentials in the differential equations for the functions $W_{n}$, $f_{n}$ become discontinuous. Since these eigenvalues cannot be used we have to devise an automatic search over the space of parameters with checks of continuity of the functions and automatic discarding of wrong results. This is accomplished by the routine P01ABF of NAG.

The search reveals that the short and respectively long radial wavelength part of the ion wave spectrum have different behavior under the modulation instability. The most affected is the long wavelength part (small $k_{x}$) with $\omega $ of the order of the ion-diamagnetic frequency. This corresponds to the dominant part of the ion spectrum, because, as can be seen in numerous numerical simulations, the dominant structures are elongated in the radial direction, with various degrees of tilting compared to the radius. It is also important to note that finding the eigenvalues in the scan of the parameter space weas possible in the condition that the eigenfunctions MATH have low number of oscillations in the radial direction (one or two). Again this corresponds to the important part of the ITG turbulence spectrum.

The variation of the parameters of the system of NSE with the physical parameters of the problem can be studied by the automatic search of eigenvalues, as described before. Below are plotted four graphs showing this dependence.
alpha1.eps
Parameter dependence of the coefficient $\alpha _{1}$.


alpha2.eps
Parameter dependence of the coefficient $\alpha _{2}$.


sigma1.eps
Parameter dependence of the coefficient $\sigma _{1}$.


sigma2.eps
Parameter dependence of the coefficient $\sigma _{2}$.