Xin Li

Bo Hu

Xieting Ling

Xuan Zeng

A Wavelet-Balance Approach for Steady-State Analysis of Nonlinear Circuits

Keywords: Nonlinear circuit, steady-state analysis, wavelet balance

In this paper, a novel wavelet-balance method is proposed for steady-state analysis of nonlinear circuits. Taking advantage of the superior computational properties of wavelets, the proposed method presents several merits compared with those conventional frequency-domain techniques. First, it has a high convergence rate ( ), where is the step length. Second, it works in time domain so that many critical problems in frequency domain, such as nonlinearity and high order harmonics, can be handled efficiently. Third, an adaptive scheme exists to automatically select proper wavelet basis functions needed at a given accuracy. Numerical experiments further prove the promising features of the proposed method in solving steady-state problems.

I. INTRODUCTION

A major difficulty in time-domain simulation of nonlinear circuits, such as power supplies, high-Q amplifiers, modulators and oscillators, etc., is that the transient response may stand for quite a long time before the steady-state is reached. This problem makes it infeasible to calculate the steady-state response by conventional transient simulation algorithms because direct integration of the circuit equations throughout the transients consumes unbearable computing time.

During the past several decades, a great number of techniques have been developed to solve the periodic steady-state problem [START_REF] Aprille | Steady-state analysis of nonlinear circuits with periodic inputs[END_REF]- [START_REF] Filicori | Nonlinear circuit analysis through periodic spline approximation[END_REF], which may be categorized into three classes: shooting methods [START_REF] Aprille | Steady-state analysis of nonlinear circuits with periodic inputs[END_REF]- [START_REF] Parkuhurst | Determining the steady-state output of nonlinear oscillatory circuits using multiple shooting[END_REF], harmonic balance methods [START_REF] Ushida | Frequency-domain analysis of nonlinear circuits driven by multi-tone signals[END_REF]- [START_REF] Yamagami | Analysis of communication circuits based on multidimensional Fourier transformation[END_REF] and sample balance methods [START_REF] Filicori | Nonlinear circuit analysis through periodic spline approximation[END_REF]. The shooting methods attempt to find a set of initial conditions satisfying the two-point boundary constraint, such that the circuit starts in periodic steady state directly. However, the shooting methods consume expensive computing time since they require to numerically integrate the system equations time after time. The harmonic balance methods assume the circuit solutions in the form of Fourier series. Moreover, they divide the circuit into a linear and a nonlinear part so that the linear subnetwork could be solved efficiently in frequency domain. Unfortunately, the harmonic balance methods need to repeatedly execute DFT and IDFT operations during the solution process, and employ a large number of harmonic components to achieve an accurate simulation result. Therefore, they also expend substantial computing time. The sample balance methods directly approximate the time-domain state-variable waveforms by suitable basis functions, such as periodic cubic splines, and use time-domain samples as problem unknowns. Nevertheless, there doesn't exist a strong theory to ensure the convergence features of these methods.

Recently, the wavelet theory has been well developed [START_REF] Meyer | Ondelettes ét opétateurs, I: Ondelettes[END_REF], [START_REF] Chui | An Introduction to Wavelets[END_REF] and widely used in many applications, such as solving partial differential equations [START_REF] Bacry | A wavelet based space-time adaptive numerical method for partial differential equations[END_REF], [START_REF] Cai | Adaptive multi-resolution collocation methods for initial boundary value problems of nonlinear PDE's[END_REF] and performing high-speed circuit simulations [START_REF] Zhou | A fast wavelet collocation method for high-speed circuit simulation[END_REF], [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF]. However, the wavelet-based method for steady-state analysis has never been explored. In this paper, we propose a novel wavelet-based algorithm, named the wavelet-balance method, to access the periodic steady-state problem. Taking advantage of the wavelets, the proposed method works more efficiently in time domain than those frequency-domain techniques, and has a solid theoretical background based on wavelet approximation theory.

The rest of this paper is organized as follows. In Section II, we first introduce the basic principle of the wavelet-balance method for steady-state analysis, then compare the proposed algorithm with conventional Fourier-based techniques [START_REF] Ushida | Frequency-domain analysis of nonlinear circuits driven by multi-tone signals[END_REF]- [START_REF] Yamagami | Analysis of communication circuits based on multidimensional Fourier transformation[END_REF] in Section III. Numerical experiments are presented in Section IV to demonstrate the computational efficiency of the proposed algorithm. Finally, we draw conclusions in Section V and outline some important background of the wavelet basis functions in the Appendix.

II. WAVELET BALANCE APPROACH FOR STEADY-STATE ANALYSIS

A. Principle of Wavelet Balanced Steady-State Analysis

Without loss of generality, we assume that a circuit is described by an ordinary differential equation of the type

dX dt = f (X; t) (1)
where [X(t) = X1(t) X2(t) 1 1 1 XN (t)] T are the N unknown state variables, and f (X; t) is a given nonlinear vector function.

The basic idea of the wavelet-balance method in this paper is to expand the unknown state variables by wavelet series, but not the Fourier series that have been used in the conventional techniques [START_REF] Ushida | Frequency-domain analysis of nonlinear circuits driven by multi-tone signals[END_REF]- [START_REF] Yamagami | Analysis of communication circuits based on multidimensional Fourier transformation[END_REF]. Assume the steady-state response period of the circuit is T , and the Sobolev space H 2 [0; L] is studied for wavelet expansion. Usually, because the period T doesn't satisfy some specific conditions needed for wavelets, for example, L 4 corresponding to the wavelet basis functions in [START_REF] Zhou | A fast wavelet collocation method for high-speed circuit simulation[END_REF], [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF], a scaling operation should be applied to map the interval [0; T ] to [0; L], i.e.,

l = K 1 t (2)
where t 2 [0; T ]; l 2 [0; L], and K is a constant which equals to L=T . Substituting (2) into (1), we obtain

K 1 dX dl = f X; l K : (3)
In order to find the steady-state solutions, the state variables X (l) are expanded by wavelets

X (l) = X 1 (l) X 2 (l) . . . XN (l) = C 11 C 12 1 1 1 C 1M C21 C22 1 1 1 C2M . . . . . . . . . . . . CN 1 CN 2 1 1 1 CN M 1 B 1 (l) B 2 (l) . . . B M (l) = C 1 B(l) (4)
where C 2 R N2M is the coefficient matrix, fBi (l); i = 1; 2; . . . ; M g, are wavelet basis functions, and M is the total number of basis functions that have been employed. The wavelet basis functions can be constructed by many means [START_REF] Chui | An Introduction to Wavelets[END_REF], but in this paper, we prefer to use the basis functions in [START_REF] Zhou | A fast wavelet collocation method for high-speed circuit simulation[END_REF], [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF] because they are proved to have a high convergence rate O(h 4 ), where h is the step length [START_REF] Zhou | A fast wavelet collocation method for high-speed circuit simulation[END_REF], [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF]. Substituting (4) into (3), we get

K 1 C 1 dB(l) dl = f CB(l); l K : (5)
Then, discretize the state (5) at some interior collocation points fl 1 ; l 2 ; . . . ; l M g [START_REF] Zhou | A fast wavelet collocation method for high-speed circuit simulation[END_REF], [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF] as shown in ( 6) at the bottom of the page.

Equations ( 1) to ( 6), are similar to those illustrated in [START_REF] Zhou | A fast wavelet collocation method for high-speed circuit simulation[END_REF] and [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF].

But for steady-state response, the state variables X(l) should further satisfy the two-point boundary constraint

X(0) = X(L): (7)
Thus, we get one more vector equation for the coefficient matrix

C C 1 B(0) = C 1 B(L): (8)
For nonautonomous circuits, the steady-state response period T is determined by input excitations, and the constant K in ( 6) is known in advance. Hence, there are M 1 N problem unknowns in coefficient matrix C to be determined. Meantime, the total number of independent algebraic equations in ( 6) and ( 8) is (M + 1) 1 N . In order to solve these equations, we define the merit function shown in ( 9) at the bottom of the page where k1k denotes the Frobenius norm. Employing the optimization algorithm, such as the Levenberg-Marquardt method in [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF], to minimize the merit function Q, we can obtain the least-squareerror solution C, as well as the steady-state response

X(l) = C 1 B(l) X(t) = C 1 B(Kt): (10)
On the other hand, in autonomous cases, there are totally M 1 N + 1 unknown variables, since the oscillation period T and consequently the constant K in (6) are undeterminated. Recall that the overall number of independent algebraic equations in ( 6) and ( 8) is (M + N ) 1 N , which is greater than or at least equals to the number of problem unknowns. Therefore, the same optimization algorithm can be used to minimize the merit function ( 9) and find the least-square-error solutions of coefficient matrix C and the unknown parameter K.

In summary, the essence of the proposed wavelet-balance method is to approximate the state variables by wavelet basis functions, but not the Fourier basis e j2k=T as in many conventional algorithms. Moreover, the wavelets are forced to keep "balance" at the two boundary points in [START_REF] Pala-Schonwalder | A discrete time approach to the steady state analysis of nonlinear autonomous circuits[END_REF], so that the steady-state solution is guaranteed when the state equation is solved.

B. Adaptive Technique

One of the main advantages of the wavelet-balance method is that there exists an adaptive scheme, which relies on the multiresolution analysis in wavelet theory [START_REF] Cai | Adaptive multi-resolution collocation methods for initial boundary value problems of nonlinear PDE's[END_REF]- [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF]. Using adaptive techniques, those wavelet basis functions, which are needed for approximating state variables, can be employed automatically. It, in turn, improves the computational efficiency significantly.

A multiresolution approximation of H 2 [0; L] is a sequence fV J ; J = . . . ; 01; 0; 1; . . .g of closed subspaces of H 2 [0; L] 1)

1 11 V 01 V 0 V 1 111 2) \ +1 J=01 V J = f0g 3) [ +1 J=01 V J = H 2 0; L] 4) VJ = VJ01 8 WJ
where the notation 8 stands for the direct sum. It is clear that the approximation accuracy depends on the wavelet space level J . The higher the space level is, the less the error will be. As pointed out in [START_REF] Cai | Adaptive multi-resolution collocation methods for initial boundary value problems of nonlinear PDE's[END_REF]- [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF], the magnitude of the wavelet coefficients in W J will indicate whether a refinement, by increasing the wavelet space level, is needed or not. For example, define the maximum relative magnitude of the wavelet coefficients in W J as

RJ = MAX C J i MAXjC i j (11)
where MAXjC J i j is the maximum magnitude of the wavelet coefficients in W J , and MAXjC i j is the maximum magnitude of all wavelet coefficients. If RJ is greater than a given error tolerance ", we increase the wavelet space level J to J 0 , where J 0 > J [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF].

More importantly, because of the compact support of wavelet bases, not all wavelet basis functions in higher wavelet spaces W J are needed in order to achieve more accuracy. In fact, only basis functions, whose positions near the singularities (which can be determined by the magnitude of the wavelet coefficients on level J ), shall be included [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF].

III. COMPARISON WITH CONVENTIONAL FOURIER-BASED TECHNIQUES

The key feature of the wavelet-balance method is the fact that wavelet basis has compact support or local support in time domain, whereas Fourier basis e j2k=T has a global support. Hence, the proposed wavelet-balance method presents following advantages.

First, the wavelet-balance method has a high-convergence rate, resulting in low computational complexity. It is demonstrated in [START_REF] Zhou | A fast wavelet collocation method for high-speed circuit simulation[END_REF], [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF] that the wavelet basis functions presented there, which are also used in this paper, have a high convergence rate O(h 4 ), where h is the step length, i.e., the distance between two adjacent collocation points.

Second, the wavelet-balance method works in the time domain, so that many critical problems in frequency domain, such as nonlinearity and high-order harmonics, can be handled efficiently. It is well known that when the circuit response is highly distorted by strong nonlinearities, the use of Fourier series with a large number of harmonics may cause considerable computation time, owing to the great number of unknown variables [START_REF] Filicori | Nonlinear circuit analysis through periodic spline approximation[END_REF].

Third, but most importantly, the adaptive scheme introduced in Section II.B can automatically select proper wavelet basis functions needed at a given accuracy. High-level wavelet basis functions are only employed near singularities. Actually, the response of nonlinear circuits will exhibit localized singularities, which can be easily captured by the wavelet approach through adaptivity. In addition, it is reported in [START_REF] Ushida | Frequency-domain analysis of nonlinear circuits driven by multi-tone signals[END_REF]- [START_REF] Filicori | Nonlinear circuit analysis through periodic spline approximation[END_REF] that the entire nonlinear network can be decomposed into a linear and a nonlinear multiport subnetwork with the same number of ports. Because the linear part can be solved in frequency domain, computation efficiency can be achieved by choosing the required number of state variables as the same number of linear subnetwork ports, no matter what the actual number of the dynamic elements is. This idea can also be applied to the proposed waveletbalance method, but the detailed implementation is beyond the scope of this paper and will not be presented here.

Q = CB(0) 0 CB(L) K 1 C dB(l1) dl 0 f CB(l 1 ); l1 K 111 K 1 C dB(lM ) dl 0 f CB(l M ); lM K 2 (9)

IV. NUMERICAL EXPERIMENTS

In this section, three circuit examples have been examined to demonstrate the effectiveness of our proposed wavelet-balance method for steady-state analysis. All of the circuit examples are simulated on a Pentium III-550 computer.

A. MOS Amplifier

Shown in Fig. 1 is a very simple MOS amplifier from which we can make a full comparison between the proposed wavelet-balance method and those conventional Fourier-based techniques. We test the circuit in Fig. 1 with a square-wave input of amplitude 61 V and frequency 1 kHz. First, we expand the state variable V OUT by Fourier series. Table I gives the relative simulation errors when different basis function numbers are employed. The relative error is defined as

Err R = [y SPICE (t) 0 y Four (t)] 2 dt [y SPICE (t)] 2 dt (12)
where y SPICE (t) represents the exact steady-state response obtained by SPICE and y Four (t) represents the simulation response by the Fourier expansion. Note that the state equation should be integrated by SPICE for quite a long time until the transients die out, so that an accurate steady-state response can be obtained. Second, we expand the state variables by wavelets and apply the wavelet-balance method to calculate the steady-state response. When the adaptive technique is not used, all basis functions in given wavelet spaces are employed for computation. Table II gives the relative simulation errors under different wavelet basis function numbers. Comparing Table II with Table I, one would find that the wavelet-balance method hasn't gained any point in this example. The convergence rate of the Fourier expansion is a little higher than that of the wavelet approach. Third, we show that the wavelet-balance method can improve its efficiency significantly after the adaptive algorithm is applied. Let the error tolerance " = 10 03 and employ the adaptive technique to se- lect proper basis functions automatically. Seven wavelet spaces (from level 0 to level 6) in all are applied by the adaptive algorithm during the course of iteration. Table III displays relative simulation error and the total number of selected wavelet basis functions as different wavelet spaces are included. Two comments can be made according to the data in Tables I-III.

1) After the adaptive algorithm is applied, the number of wavelet bases which are used at a given space level is reduced since not all basis functions in that wavelet space are needed for computation. On the other hand, the simulation error for given space level J doesn't change much as shown in Tables II andIII. This demonstrates that the adaptive algorithm has the potential to pick up the most important basis functions and those ones neglected by the adaptive scheme don't contribute a lot to the simulation accuracy. 2) Taking advantage of the property of compact support, the adaptive algorithm can automatically select proper wavelet basis functions. However, Fourier bases have global support and it is infeasible to realize such a selection in time domain as wavelets. Therefore, as displayed in Tables I and III, the wavelet expansion exhibits higher convergence rate than the Fourier approach after the adaptive scheme is employed. Finally, Fig. 2 depicts the relation between simulation accuracy and computational time for the above three different approaches. Note that under small basis function numbers, the wavelet-balance method is not so efficient as the Fourier approach because almost all low-level wavelet basis functions are employed in order to capture those globalized components of the circuit response. However, as the wavelet space-level increases, the adaptive scheme begins to take effect and high level wavelet bases are only used in those regions where the circuit response exhibits localized singularities. Therefore, under large basis function numbers, the wavelet-balance method begins to catch up with the Fourier approach in accuracy and become more efficient.

B. Power Supply

The power-supply circuit shown in Fig. 3 is one of the several examples that has been used in the past to test various steady-state analysis algorithms for nonautonomous circuits [START_REF] Aprille | Steady-state analysis of nonlinear circuits with periodic inputs[END_REF], [START_REF] Skelboe | Computation of the periodic steady-state response of nonlinear circuits by extrapolaton methods[END_REF], [START_REF] Parkuhurst | Determining the steady-state output of nonlinear oscillatory circuits using multiple shooting[END_REF], [START_REF] Frey | An integral equation approach to the periodic steady-state problem in nonlinear circuits[END_REF]. The input excitation in Fig. 3 a standard sinusoidal signal of amplitude 610 V and frequency 60 Hz. First, we expand the state variables V C1 ; V C2 ; V C3 and I L4 by Fourier series. Table IV displays the simulation errors and computation times when different basis function numbers are used. Second, we expand the state variables by wavelets and apply the adaptive scheme to automatically choose necessary wavelet basis functions with error tolerance " = 10 03 . Table V gives the simulation results obtained by the wavelet-balance method. The effect of the adaptive algorithm is well illustrated by Tables IV and V that the wavelet expansion overtakes the Fourier technique in efficiency as the basis function number increases.

Figs. 4 and 5 depict the steady-state response V OUT and I L4 ob- tained by the wavelet-balance method at space level J = 3. As indicated by the distribution of collocation points 1 , the wavelet basis functions which are adaptively employed for representing IL4 are much more than those for V OUT , since the waveform of I L4 is more singular than that of V OUT . In addition, most wavelet basis functions in Fig. 5 concentrate in interval [0, 0.005] where the response of IL4 exhibits a fast transition. The above comparison demonstrates that the waveforms with different singularities are treated by the adaptive scheme differently. High-level wavelet bases are only used near sharp transitions, which is exactly the reason why the wavelet-balance method is more efficient than the Fourier approach in this example.

C. Van der Pol Oscillator

Because of the high nonlinearity, the Van der Pol oscillator is served as a benchmark circuit for various steady-state algorithms for autonomous circuits, such as those in [START_REF] Parkuhurst | Determining the steady-state output of nonlinear oscillatory circuits using multiple shooting[END_REF], [START_REF] Pala-Schonwalder | A discrete time approach to the steady state analysis of nonlinear autonomous circuits[END_REF], [START_REF] Frey | A class of relaxation algorithms for finding the periodic steady-state solution in nonlinear system[END_REF]. The circuit schematic is quite simple as illustrated in Fig. 6. Tables VI and VII display the simulation results obtained by the Fourier expansion and 1 Each wavelet basis function is combined with a collocation point, and the wavelet basis function is localized in the neighborhood of its corresponding collocation point. Therefore, the distribution of collocation points reflects the distribution of wavelet basis functions. the wavelet-balance method (with adaptive scheme at error tolerance " = 10 02 ) respectively. Again, the wavelet-balance method reveals its effectiveness in both speed and accuracy when solving the steady-state problems. Fig. 7 depicts the output waveform V OUT simulated by the waveletbalance method at space level J = 4. Note that the adaptive scheme only automatically employ high level basis functions in interval [START_REF] Skelboe | Computation of the periodic steady-state response of nonlinear circuits by extrapolaton methods[END_REF][START_REF] Ushida | Frequency-domain analysis of nonlinear circuits driven by multi-tone signals[END_REF] and [START_REF] Frey | An integral equation approach to the periodic steady-state problem in nonlinear circuits[END_REF][START_REF] Yamagami | Analysis of communication circuits based on multidimensional Fourier transformation[END_REF], where fast transitions take place.

V. CONCLUSION REMARKS

We present in this brief, a novel wavelet-balance method for steadystate analysis of nonlinear circuits. As a counterpart of those Fourierbased techniques, the wavelet-balance method works in time domain and has a solid theoretical background based on the wavelet approximation theory. Taking advantage of the superior computational properties of wavelets, the proposed method can automatically select proper basis functions for computation so that it presents greater efficiency than the conventional Fourier-based techniques in many practical problems. From this point of view, the wavelet-balance method exploits a new approach to access the steady-state problem besides those frequency domain methods that have been used for a long time.

K 1 C

 1 1 dB(l1) dl dB(l2) dl 111 dB(lM ) dl = f CB(l 1 ); l1 K f CB(l 2 ); l2 K 11 1 f CB(l M ); lM K (6)

Fig. 1 .

 1 Fig. 1. MOS amplifier where = 1 k and = 0 1 F.

Fig. 3 .

 3 Fig. 3. Power supply where

Fig. 4 .

 4 Fig. 4. Steady-state response of in power supply.

Fig. 5 .

 5 Fig. 5. Steady-state response of in power supply.

Fig. 6 .

 6 Fig. 6. Van der Pol oscillator where = 1 F, = 1 H and ( ) = 5 ( 3).

Fig. 7 .

 7 Fig. 7. Steady-state response of in Van der Pol oscillator.

TABLE I SIMULATION

 I RESULTS OF THE FOURIER EXPANSION FOR MOS AMPLIFIER

TABLE II SIMULATION

 II RESULTS OF THE WAVELET EXPANSION (WITHOUT ADAPTIVE SCHEME) FOR MOS AMPLIFIERTABLE III SIMULATION RESULTS OF THE WAVELET EXPANSION (WITH ADAPTIVE SCHEME)

FOR MOS AMPLIFIER Fig. 2. Simulation results for MOS amplifier.

TABLE IV SIMULATION

 IV RESULTS OF THE FOURIER EXPANSION FOR POWER SUPPLYTABLE V SIMULATION RESULTS OF THE WAVELET EXPANSION (WITH ADAPTIVE SCHEME) FOR POWER SUPPLY

TABLE VI SIMULATION

 VI RESULTS OF THE FOURIER EXPANSION FOR VAN DER POL OSCILLATOR TABLE VII SIMULATION RESULTS OF THE WAVELET EXPANSION (WITH ADAPTIVE SCHEME) FOR VAN DER POL OSCILLATOR

ACKNOWLEDGMENT

The authors wish to thank Prof. D. Zhou, University of Texas, Dallas for the helpful discussions, and all the anonymous reviewers of this paper for their valuable suggestions.

APPENDIX

In this Appendix, we briefly present some important background of the wavelet basis functions used in this paper.

In Sobolev space H 2 [0;L], we first introduce the following function subspaces:

V 01 = span 1(t); 2(t); 1(L 0 t); 2(L 0 t); '0;01(t) ' 0;0 (t);' 0;1 (t);.. .;' 0;L04 (t);' 0;L03 (L 0 t) W J = spanf9 J;K (t);01 K 2 J 0 2g; J 0 V J = V 01 + W J ; J 0

where the definitions of functions 1 (t); 2 (t);' 0;01 (t);f' 0;K (t); 0 K L 0 4g;' 0;L03 (L 0 t) and f9 J;K (t);01 K 2 J 0 2;J 0g are given below. spanff1;f2; .. .;fNg is a function set formed by all linear combinations of the functions f 1 ; f 2 ; .. .;f N .

In (A.1), the function 1 (t); 2 (t) are used to handle nonhomogeneity of the boundary data

where

The boundary scaling functions are The boundary wavelet functions are and the interior wavelet functions are J;K (t) = (2 J 0 K); 1 K 2 J 0 4

(A.9) where (t) = 0 3 7 '(2t) + 12 7 '(2t 0 1) 0 3 7 '(2t 0 2):

(A.10)

The total number of basis functions in subspace VJ is N = 2 J +1 L+3.

On the other hand, the following collocation points are chosen for V 01 and W J ; J 0, respectively t = 0; 1 2 ; 1; 2; . ..;L 0 1;L 0 1 2 ; L;

for V 01 t = 1 2 J +2 ; 1:5 2 J ; 2:5 2 J ; ...; 2 J 0 1:5 2 J ; L 0 1 2 J +2 ; for W J (A.11)

The above content only gives a simple outline of the wavelet basis functions and collocation points used in this paper, and further details can be found in [START_REF] Cai | Adaptive multi-resolution collocation methods for initial boundary value problems of nonlinear PDE's[END_REF]- [START_REF] Zhou | An adaptive wavelet method for nonlinear circuit simulation[END_REF].


Professor. W J H K Wang

Chen

Investigation of the Vibration of a Blade With Friction Damper by HBM

The friction damper has been widely used to reduce the resonant vibration of blades. The most commonly used methods for studying the dynamic behavior of a blade with a friction damper are direct integration methods. Although the harmonic balance method (HBM) is a well-known method for studying nonlinear vibration problems, generally only a one-term approximation has been proposed to study the nonlinear vibration of a frictionally damped blade. In this work, a HMB procedure with a multiterm approximation is proposed. The results show that the steady-state response and other related behavior of a frictionally damped blade can be predicted accurately and quickly by an HBM with a multiterm approximation.

Introduction

The friction damper has been widely used to reduce the resonant vibration of blades. In the past, many works have been done to investigate the dynamic behavior of a blade with a friction damper either by the macroslip approach [START_REF] Griffin | Friction Damping of Resonant Stress in Gas Turbine Engine Airfoil[END_REF][START_REF] Srinivasan | Dry Friction Damping Mechanisms in Engine Blades[END_REF][START_REF] Wang | Design of Shroud Interface-Angle to Minimize the Forced Vibration of Blades[END_REF] or by the microslip approach (Menq et al., 1986a(Menq et al., , 1986b)). Generally, the steady-state nonlinear forced vibration of a blade was investigated by direct time-step integration with iterative process. It is very time consuming because the decay of the transient vibration generally is very slow. Although the harmonic balance method (HBM) is a well-known method to study the nonlinear vibration problems, usually only the HBM with one-term harmonic was used to investigate the harmonic forced vibration of a blade with a friction damper [START_REF] Menq | A Comparison of Transient and Steady State Finite Element of Forced Response of a Frictionally Damped Beam[END_REF]. This is due to the difficulty that the nonlinear friction force is dependent on the stepwise stick-slip motion, which is not known in advance. As a consequence, the number of simultaneous algebraic equations found from the HBM is less than the number of unknowns if the response is assumed to consist of many harmonic terms. The slowly varying parameter method has also been used to study the harmonic response of a blade with a friction damper [START_REF] Sinha | Effects of Static Friction on the Forced Response of Frictionally Damped Turbine Blades[END_REF]. However, only one term of harmonic has been considered in that work. A modified HBM with a multiharmonic approximation has been proposed by [START_REF] Pierre | Multi-harmonic Analysis of Dry Friction Damped Systems Using an Incremental Harmonic Balance Method[END_REF] to investigate a dry friction damped system. The method is also known as incremental harmonic balance (IHM). However, the IHB method is no more efficient than direct numerical integration when the number of harmonics is increased. Besides the IHB, an alternating frequency/time (AFT) method [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response on Nonlinear Dynamic Systems[END_REF] the study of a dry friction problem. However, the method has to switch the calculation alternately in frequency and time domains. Besides, the aliasing and leakage problems also must be handled carefully.

The traditional HBM due to its simplicity has been widely accepted in the past. However, as mentioned, it is difficult to apply the traditional HBM with multiple harmonics to investigate the vibration of a dry friction damped blade. In this work, a method was developed to overcome the difficulty when the HBM with many harmonic terms was used. Some extra constraint equations were found to supplement the algebraic equations found from the traditional harmonic balance method. The accuracy and the computational time from the multiterm HBM were compared with the results from the direct time integration.

Analysis

As mentioned, two different friction models, i.e., macroslip and microslip approaches, have been proposed in the past. Generally speaking, the microslip approach can provide more accurate results only when the normal preload of the damper is high, and at the expense of higher computational effort, as discussed by Menq et al. (1986a) and [START_REF] Wang | The Influence of a Variable Friction Coefficient on the Dynamic Behavior of a Blade With a Friction Damper[END_REF]. Therefore, the macroslip approach was adopted in this work.

A typical single blade with a massless, flexible, blade-toground friction damper is shown in Fig. 1. If the blade is discretized by the finite element method, then the equation of motion of the blade can be written as: single degree of freedom system, as shown in Fig. 2. The equation of motion of the single-mode model can be written as

[M\{X}+[C]{X] + [K]{X) = {F} + {F n ) ( 1
mx + ex + kx =f-f n (2)
where / is the external excitation force and /" is the friction force due to the friction damper. N, /x, and k d in Fig. 2 represent the normal preload, the friction coefficient of the friction damper, and the stiffness of the damper in the direction of relative motion, respectively. According to the macroslip approach, the friction force, f", can be expressed as

( + k d {x-y) when k d \x-y\<y.N ]^ +liNsign(y) when k d lx-y\>fiN
where y is the displacement of the friction damper, as shown in Fig. 2. If the external excitation force is a simple harmonic force f=fo x cos wt, the steady-state solution of Eq. ( 2) can then be found by the harmonic balance method (HBM). In the following section, the one-term, two-term, and multi-term approximations will be discussed separately.

One-Term Approximation. According to the HBM the steady-state solution of Eq. ( 2) can be expressed as a Fourier expansion,

x= ^TJ [Ajcos (iwt) +-B,-sin (iwt)} i'=i (4)
In this section, only a one-term approximation will be used, i.e.,

x=A t cos wt + B t sin wt = A l cosd + B l sm9_ = Rcosd (5)
where 9 = wt and 0 = w/ -</ >. The angle <$> is the phase difference between the external excitation force and the response. According to Eq. ( 5) the angle 9 should be equal to zero when the response is maximum, R. A typical relationship between the response and the friction force is shown in Fig. 3. The damper is in stick and slip conditions during the AB and BC sections, respectively. Then the velocity of the damper, y (indicated in Fig. 2), should be: fO whenO<0<0* (6) y = x when 0*<0S7r From Eq. ( 6), one has

(x + c 2 whenO<0<0* when 0*<0<ir (7)
The two constants c x and c 2 can be determined by the analysis in what follows. During the stick section, the friction force is equal to zero when x is equal to>>. According to Fig. 3,/" = 0 (the middle point of line AB), when x is equal to R -fxN/k d .

In other words, one finds that c\ = R -nN/k d . During the slip section, the difference between x and y is held constant, i.e., c 2 which is the horizontal distance between point B and the middle point of line AB in Fig. 3. Then one can find that c 2 = fiN/k d . From Eq. ( 3), the friction force can then be expressed as:

f n = k d (x-y) = k d (x-R)+nN, whenO<0<0* (8«) (8b) f n = ltNsign(y) = liNsign (x) = -nN, when0*<0<ir
Equation (8a) should be equal to Eq. (&b) when 6 = 9*, then one has

6* = cos~x{\ -2fiN/(k d R)] (9)
Because the response is assumed to be a simple harmonic form, the response cycle can be divided into two symmetric half cycles, i.e., the half cycle ABC is equal to the other half cycle CD A in Fig. 3. So, one can find that /"(» + *)=-/"(»), O<0<TT

With Eqs. ( 8) and ( 10), the friction force can be expanded by Fourier series, and only the first harmonic term is kept, /"(0)=F c cos0 + F s sin0 (11)

with

F c = ( -k d R/Tr)[d* -0.5sm(2d*)} F s = (-VV/TT ) [ 1 -pN/ (k d R) ]
Substituting Eqs. ( 5) and ( 11) into Eq. ( 2) and separate the coefficients of cos0 and sin0, one has: 12c) and (12c/) are nonlinear equations with two unknowns R and <j>, which can be solved by the traditional Newton-Raphson method. In practice, only the unknown R needs to be found by iteration; the angle <j> can be obtained directly from Eq. (12c?) providing the R is known. The response is then known from Eq. ( 5). The one term approximation is not the target of the present work; however, it may provide valuable information for multiterm approximation. Two-Term Approximation. From the result of the direct time integration, it is known that the steady-state response of the system in Fig. 2 contains only the odd harmonic components of the external excitation force. Therefore, if a two-term approximation is assumed, according to Eq. ( 4), the response can be written as x=Ai cos w t + B x sin wt + A 3 cos 3 wt + B 3 sin 3 wt (13) = y4icos0 + j 8isin0+yl3cos30 + 53sin30 Let the maximum value of x be denoted by A m and the corresponding 0 be denoted by 0 O , as shown in Fig. 4. Point B in Fig. 4 is the slip impending point, and the value of 0 at this point is indicated by 0*. In a manner similar to the one-term approximation, the friction force during the stick portion (line AB) can be expressed as

or \(k- (k-mw 2 )R+F c =f 0 cos<f> cwR -F s =/ 0 sin </ > mw 2 )R + F C ] 2 + [cwR -FJ 2 =/ 0 2 <t> = sin">[(cwR-F s )/f 0 ] (12a) (12b) (12c) (12rf) Equations (
f" = k d (x-A,") + fiN = kalA 1 (cos6-cosd 0 ) + Bi(sm6 -smd 0 ) + A 3 (cos36_-cos38o) + 5 3 (sin30-sin30o)] + /JV, when 0 O <0<0* (14)
During the slip portion, the friction force is /"= -iiN, when0*<0<7r + 0 o (15) So far, 0 O and 0* are unknown. This is the main difference between the one-term and the multiterm approximations. How to determine these unknowns will be discussed later. With Eqs. ( 14) and ( 15), the friction force can be expanded as a Fourier series; however, only the first two terms are kept,

f" = F cl cos 6+F sl sin 6+ F ci cos 36+ F s3 sin36 (16)
The coefficients -F cl , F sl , F c3 , and F s3 can be found in the appendix. By substituting Eqs. ( 13) and ( 16) into Eq. ( 2), and setting the coefficients of cos0, sin0, cos 30, and sin 30 equal to zero, one has If 0 O and 0* are known, then Eqs. ( 17)-( 20) are linear equations with A u A 3 , B u B 3 as unknowns, which can be solved easily.

(k-mw^Ai + wcBi-fo + Fc^O . ( 17
) (k-9mw 2 )A 3 + 3wcB 3 + F c3 = 0 (18) {k-mw 2 )B l -wcA l +F sl = Q (19)
The values of 6 0 and 0* can be determined by Eqs. ( 21)-( 23), which are nonlinear. In summary, Eqs. ( 17)-( 23) can be used to determined the unknowns A it A 3 , B u B 3 , 0 O , 0*. The procedure is summarized as follows:

1 Give proper initial guess values of 0 O and 0* and then F cl , F c3 , F sl , and F s3 can be found (see appendix).

2 The unknowns A it A 3 , B u B 3 can be determined from the linear equations Eqs. ( 17)-( 20) by direct matrix inversion.

.3 With A U A 3 , B\, and£ 3 , the condition of Eq. ( 22) is then checked:

(a) If Eq. ( 22) is satisfied, then go to step (4).

(b) If Eq. ( 22) is not satisfied, then the initial guess values of 0 O and 0* are modified as: new 0 o =old 0o+7r, new 0* = old 6* + w. After the modification, then go back to step (1).

4 With A u A 3 , B u and B 3 , Eqs. ( 21) and ( 23) are then used to find the exact values of 0 O and 0* by the Newton-Raphson iteration method. The flow chart of the procedure is given in Fig. 5.

The modification in step (3) can be understood easily be- It is known that the most important step in using the Newton-Raphson iteration method is the guess of the initial values. How to guess proper initial values of 6_ 0 and 0* will be discussed in the next section.

The above derivation and solution procedure can easily be extended to the multiterm approximation. If an «-term approximation is used, then one can find 2/z linear equations like Eqs. ( 17)-( 20) and two nonlinear equations like Eqs. ( 21) and ( 23).

Results and Discussion

In this section, three main topics will be discussed, i.e., the accuracy of the HBM, the computation time, and the effect of the initial guess values. Before the discussion, the system parameters of Fig. 2 are given here: m = 1.24, k = 17890, c= 3, and k d = 3500.

Accuracy of the HBM. As mentioned, generally only oneterm HBM was used to investigate the vibration of a blade with a friction damper. In some cases, for instance, in fully stuck or fully slipping conditions, the one-term approximation is accurate enough to predict the response. However, the oneterm approximation generally may cause significant error.

Here, an example is given to show this situation. Figure 6 shows a response found by the direct time integration, one- term HBM, two-term HBM and three-term HBM, respectively. The vertical axis is the normalized response A/A m , where A m is the maximum amplitude found by the direct time integration. In this case the excitation frequency w is set to be 1/3 w 0 , where w 0 is the natural frequency of the system of Fig. 2 an optimal normal load TV, which results in minimum response, can be found. The normalized normal load ixN/f 0 in Fig. 6 is set to the optimal value pN/f 0 = 0.145. One can find that the one-term approximation results in significant error. In this case, the three-term approximation can accurately predict the response, and its result is difficult to distinguish from the result of the direct integration in the figure. The corresponding frictional damping loop is shown in Fig. 7. One can find that the frictional damping cannot be simulated accurately by the oneterm approximation. The above result indicates that only the multiterm HBM can accurately predict the response of a blade with a friction damper. The result also demonstrates that the multiterm HBM proposed in this work is feasible. Theoretically, the HBM with more terms will result in a more accurate solution. In practice, according to our experience, it needs no more than three terms. The question of the computational time will be discussed in the next section.

Computer Time. No matter how many terms, say n terms, are kept in the HBM, only two nonlinear equations (i.e., Eqs. ( 21) and ( 23)) and 2n simultaneous linear equations (i.e., Eqs. ( 17)-( 20) are needed to find all the unknowns. As shown in Fig. 5, the two nonlinear equations are solved by the iteration method, and the 2n linear equations must be solved once for each iteration. In other words, if m iterations are needed until convergence is reached, then the In linear equations must be solved m times. Therefore, how to reduce the number of iteration is very important. It is known that the initial guess values are very important for the iteration process. An improper initial guess value may drastically increase the number of iterations or cause divergence. Thus, to discuss the computational time, the effect of the number of terms should be considered together with the effect of the initial guess values. First, the effect of the initial guess values will be discussed.

The nonlinearity of the vibratory blade comes only from that nonlinear friction force. It is known that the friction force can be expanded by Fourier series with the fundamental frequency equal to the external excitation frequency w. In other words, the vibration of the blade is mainly caused by a harmonic force with a frequency equal to w. Thus, the solution obtained from the one-term approximation should be a good initial guess value for the multiterm approximation. For the multiterm approximation, one should give the initial guess values of 8 0 and 8*. In this work, it is suggested that the 6 0 and 6* obtained from the one-term approximation can be used as the initial guess values of 8 0 and 8* for the multiterm approximation. According to Eq. ( 5), one has the relations: 8 0 = 8 0 + <p and 8* = 8* + <j). As pointed out in the one-term approximation, do is zero, while 8* and <j> can be obtained from Eq. ( 9) and Eq. ( 12), respectively. To understand the effect of the initial guess values on the computer time, different guess values are assumed and the necessary computer time is recorded. The result is shown in Fig. 8. The horizontal axis "shift of initial value" means the amount the initial guess values of 6 0 and 8* are shifted from the suggested initial values, i.e., 8 0 = 4> and 8*=8* + 4>. For instance, "-0.2" indicates that initial guess values are given as: £ 0 = (1 -0.2)<^>, 8* = (1 -O.2)(0* + 0). One can find that if the shift of the initial value is between -0.5 and 0.3, one can find a converged solution corresponding to the solution found from the direct integration. If the initial guess values are beyond this range, one may find another converged solution or the solution is diverged. Here "another converged solution" means a solution different from the solution found from the direct time integration. So far, it is not certain that this "another solution" exist in practice. The result of Fig. 8 indicates that the 8* and 4> found from one-term approximation are the proper initial guess values of 8 0 and 8* for the multiterm approximation.

If the initial guess values are selected properly, the ratios of the computer time are 1:2.4:5.5:61 for one-term, two-term, three-term approximations and the direct time integration, respectively. Note that the computer time for the two-term and three-term approximations includes the portion to find the proper initial guess values, i.e., 8* and 0. As mentioned, the HBM requires no more than three terms in order to obtain an accurate result. That is to say the HBM needs a computer time about one order of magnitude smaller than that of the direct integration. The above result is concluded not only from the above example, but also from many other calculated examples.

The advantage of the proposed HBM is especially significant when the response must be recalculated many times, for instance, to find the optimal normal load or the frequency response curve. It is well known that there exists an optimal normal load TV with which the blade will experience minimum forced vibration. The optimal normal load is an important design parameter of the friction damper. Generally, the optimal normal load cannot be found directly except when the response is a pure simple harmonic function. The optimal normal load generally is found by observing the response for many different normal load values. In other words, a curve indicating the responses versus the normal load TV should be created, and the minimum response corresponds to the optimal normal load. For each value of TV, a new calculation must be executed. It is very time consuming if the direct integration method is used. For the proposed method, the answer can be found very quickly, because the solution for normal load TV can be used as a very good initial guess value for finding the solution at TV+ ATV; ATVis a small increment. Because the initial guess values are very close to the exact values, the iteration converges very quickly. The frequency, response curve, i.e., the response versus the excitation frequency, can also be found quickly by the proposed method. The response at a certain excitation frequency w can be used as a very good initial guess value for finding the response at frequency w + Aw.

In summary, there are two important features of the proposed method: (1) There are only two nonlinear equations no matter how many terms are kept in the HBM; (2) a proper initial guess value for the iteration generally can be found easily. Due to these features, the proposed method is very efficient for the steady-state response and related analyses.

Conclusion

The dynamic behavior of a blade with a friction damper has been intensively investigated in the past. Although the harmonic balance method (HBM) is a well-known method, usually only the HBM with a one-term (or one harmonic component) approximation was used to investigate the harmonic forced vibration of a blade with a friction damper. Generally, the HBM with an n-term approximation will result in 2« simultaneous nonlinear algebraic equations in terms of 2« unknown Fourier coefficients for a system with one degree of freedom. However, this is not the case in the proposed method. The main step of the proposed procedure is to find the impending time points (i.e., 0 O and 0*) of stick and slip of the friction damper. Because the 0 O and 0* are found from two nonlinear algebraic equations, the problem of the initial guess values of the do and 0* has also been discussed in this work. The results show that the HBM needs no more than three terms in order to obtain an accurate result. The main advantage of the HBM is the computational efficiency; it needs only about one-tenth of the computer time in comparison with the direct integration. Due to the drastic reduction of the computer time, it becomes possible to investigate many important characteristics of a frictionally damped blade in a very short time.

  can also accommodate multiple harmonics for Contributed by the International Gas Turbine Institute and presented at the 37th International Gas Turbine and Aeroengine Congress and Exposition, Cologne, Germany, June 1-4, 1992. Manuscript received by the International Gas Turbine Institute January 16, 1992. Paper No. 92-GT-8. Associate Technical Editor: L. S. Langston.

  Fig. 1 Schematic of a single blade with a friction damper

Fig. 4 -

 4 Fig. 4 Frictional damping loop, multiterm approximation(k-9mw 2 )B 3 -3wcA 3 +F s3 = 0 (20)Because the F cU F si , F c3 , and F s3 are functions of 0 O and 6* which are unknown, the four equations Eqs. (17)-(20) have six unknowns, i.e., A u B iy A 3 , B 3 , 0 O and 0*. So, one must find other equations to supplement Eqs. (17)-(20). The response x is maximum (with value A,") when 0 = 0 O . Thus, the following equations should be satisfied:dx -= -A i sin 0 O + -Si cos 0 O at

Fig. 5

 5 Fig. 5 Flow chart of the proposed iteration procedure

  Fig. 6 Steady-state forced response, w/w0 = 0.334, pN/f0 = 0.145

  Fig.7Frictional damping loop; the legend is the same as Fig.6

294/Vol. 115, APRIL 1993Transactions of the ASME

298/Vol. 115, APRIL 1993Transactions of the ASME

APPENDIX

The coefficient of the Fourier series in Eq. ( 16) is derived in this appendix. + -Bi(sin0-sin0o) + 5 3 (sin30-siri30o)]cos0rf0

IxNcos 36d6


X Y Leng

D Krasnov

Y Kolesnikov
email: yuri.kolesnikov@tu-ilmenau.de

B W Li

The distinction of instability and turbulence between the Taylor-Couette flow and electrically driven flow in annular channel

Keywords: Taylor-Couette flow, electrically driven flow, turbulence

. The approximation of low magnetic Reynolds number, corresponding to liquid metal flows, is used. In the case of TCF, the shear flow is driven by inner cylinder, and the EDF is driven by the mean Lorentz force, which is a product of applied field and electric currents injected in the radial direction. The latter mimics the case of pressure driven flows. Therefore, despite similar geometry, the mechanisms of generation of shear and turbulence are different, as well as potential routes to the suppression of eddies by magnetic field. The resemblances and differences of turbulent behaviors of both flows are revealed in present study.

Introduction

Taylor-Couette flow (TCF, hereafter) is the typical shear flow driven by rotating inner cylinder; analogous to the Dean flow, electrically driven flow (EDF, hereafter) is the typical Poiseuille flow in which the pressure gradient is replaced by the Lorentz force. For the both flows, the electrical conducting fluid is confined between two concentric cylinders with the same geometry, e.g. aspect ratio and radius ratio. Apparently, the mechanisms of these two flows are distinct. However, unlike the plane Couette flow and the plane Poiseuille flow, due to the presence of the centrifugal instability in concentric cylinders, there will be more different and interesting turbulent behaviors. Meanwhile, comparing these two flows in the same geometry, there also will be some similar and dissimilar behaviors. Via the comparative study of Taylor Couette flow and electrically conducting flow, we hope to reveal the differences of the turbulent behaviors. We consider the flow of the incompressible viscous electrically conducting fluid with constant kinematic viscosity  , electrical conductivity  , magnetic diffusivity  , magnetic permeability 0  , and density  , contained in the annulus between two concentric cylinders of radii 1 r and 2 r . The flow is subject to an axial uniform magnetic field 

 z B B e ,
where z e is the unit vector in axial direction.

   ).
Thereafter, the governing equations written in dimensionless form are,

    2 2 1 ()                       u Ha u u p u u e e t Re Re (1) 0    u (2) 2 ()       ue (3) ()        j u e ( 4
)
where t is the dimensionless time,  u ( u , v , w ) are the dimensionless velocities components in the azimuthal, axial and radial direction, respectively, p is the dimensionless pressure,

 e ( 0   e , 1  z e , 0  r e
) are the components of the magnetic field in three directions respectively, for cylindrical coordinates (  , z , r ). In the present study, we use periodical boundary condition in axial direction, and the dimensionless domain length is 2

 .The radius ratios are 0.8   in the both flows. In TCF, the potential of both cylinders are 0; and in EDF, the potential of inner and outer cylinder are 0 and 1 respectively. For velocity, no-slip boundary condition is applied for both cylinders.

Numerical method

The equation system (1-4) is solved numerically using the finite difference scheme via direct numerical simulation. The scheme is of the second order of approximation and is based on the spatial discretization, which is nearly fully conservative in regard of mass, momentum, kinetic energy, and electric charge conservation [3,4]. The explicit Adams-Bashforth/Backward-Differentiation method of the second order is employed for the time discretization. The velocity and pressure are by the projection method. At every time step, the two Poisson equations-the projection method equation for pressure and the equation for potential (3) are solved using FFT in the azimuthal direction and cyclic reduction direct solver(a part of the software package Fishpack) in the meridional plane. The OpenMP interface for shared memory multiprocessors is used for the parallelized algorithm. More details of the numerical method are described in [2]. It is possible for us to compare the two results. The velocity fluctuations are defined according to the average in time, in axial and azimuthal direction [5],

Results

As the definitions of

   zt u u u (5)
where  u is the velocity fluctuation due to the turbulence and mean Taylor vortex motion(or mean azimuthal vortex motion in EDF). The velocity fluctuations are defined according to the average in time and azimuthal direction,

   t u u u (6)
where  u is the velocity fluctuation due to the turbulence. We compare the root mean square velocity fluctuations in three directions and the kinetic energy of velocity fluctuations in two flows. Kinetic energy of velocity fluctuations of the turbulence and mean Taylor vortex motion can be written as

Conclusions

In the present study, with the radius ratio 0.8   and 1800

 mean Re , we compare the root mean square velocity fluctuations and kinetic energy of fluctuations in Taylor-Couette flow and electrically driven flow. The contributions from the mean TV motion (or vortex flow in EDF) play important roles in these two turbulent flows, and the contributions of mean flow in TCF are much larger than that in EDF. The mean motion in EDF has less effect on the turbulence. For the both flows, the velocity fluctuations of the azimuthal and axial velocity near the wall are stronger than those in the middle region. For the EDF the turbulent fluctuations near the outer cylinder are larger than that near the inner cylinders; for the TCF, turbulent fluctuations near the cylinders are same.

  Re are not same. So, we use the similar mean Re based on the mean azimuthal velocity as the standard for the comparison. The cases that we consider here are the TCF for 4000

Fig. 1 :Fig. 2 :Fig. 3 :

 123 Fig. 1: r. m. s. velocity fluctuations in TCF.

Acknowledgment

We are grateful to Christian Karcher, Thomas Boeck and Andre Thess, for their interesting discussions. Financial support was provided by the Deutsche Forschungsgemeinschaft (DFG) within the RTG Lorentz Force Velocimetry and Eddy Current Testing. X.Y. Leng was supported by the China Scholarship Council and the NSFC (Grant Nos. 10772044 & 51176026). Computer resources were provided by the computing center of TU Ilmenau.


A Pothérat
email: alban.potherat@coventry.ac.uk

Quasi-static magnetohydrodynamic turbulence between two and three dimensions

Keywords: MHD, turbulence dimensionality, two-dimensional turbulence, turbulence promoters

Recent concepts of quasi-static MHD turbulence between plane walls are presented. It is shown that the dimensionality of this type of flow is governed by the ratio of the diffusion length associated to the Lorentz force to the channel width. Depending on turbulence dimensionality, three different dissipation mechanisms are activated that correspond to three different scalings for the intensity of turbulent fluctuations. In all three regimes, the relative turbulent intensity is found to increase with the applied magnetic field, in apparent contrast to commonly accepted ideas.

Introduction

Turbulence under an externally imposed, static magnetic field has a well-known tendency to two-dimensionality [START_REF] Moffatt | On the suppression of turbulence by a uniform magnetic field[END_REF]. This effect is driven by the diffusive nature of the Lorentz force at low magnetic Reynolds number [START_REF] Sommeria | Why, how and when mhd turbulence becomes two-dimensional[END_REF] which smoothes out velocity gradients along the magnetic field. The final dimensionality of the flow, is however strongly determined by the nature of the boundaries of the fluid domain. In the generic case of a channel perpendicular to the field and bounded by two walls, Hartmann boundary layers that develop along the walls preclude a fully quasi-two dimensional state, so the flow can be at best twodimensional in its bulk, or quasi-two-dimensional [START_REF] Sommeria | Why, how and when mhd turbulence becomes two-dimensional[END_REF][START_REF] Pothérat | Mhd turbulence at low rm: the role of boundaries[END_REF]. With non dissipative boundaries (periodic or slip-free boundary conditions), the transition between strictly two-dimensional and three-dimensional states is mainly governed by the stability properties of large two-dimensional structures [START_REF] Boeck | Large-scale intermittency of liquid-metal channel flow in a magnetic field[END_REF]. When walls are present, by contrast, the flow states span a continuous spectrum of states involving different types of three-dimensional effects (presence of transversal velocity, weak velocity gradients preserving topological equivalent between planes perpendicular to the magnetic field, full three-dimensionality) [START_REF] Klein | Appearance of three-dimensionality in wall bounded mhd flows[END_REF][START_REF] Pothérat | Why, how and when mhd turbulence at low rm become three-dimensional[END_REF][START_REF] Baker | Three-dimensional stucture of electrically driven vortices[END_REF]. Considering a generic channel configuration, we show that these states of dimensionality are solely determined by the ratio l z /h of the diffusion length associated to the Lorentz force to the channel width. We show that three different states can be distinguished, each of which characterised by a how the intensity of turbulent fluctuations scales with the externally applied forcing. These scalings imply that in the process of making the flow quasi-two dimensional, the effect of the magnetic field is to lower Joule dissipation to the point where it actually increases the intensity of turbulence. These theoretical concepts are verified experimentally on the FLOWCUBE experimental platform where turbulence is driven in a cubic container by injecting electric current at one of the Hartmann walls of the vessel. In FLOWCUBE, the intensity of the forcing is measured directly by the quantity of current that is injected in the flow.

Dimensionality and turbulence intensity

Consider a horizontal channel of width h, filled with liquid metal (density ρ, viscosity ν, electric conductivity σ) and pervaded by a homogeneous static magnetic field Be z . Both the tendency to twodimensionality and the sources of three-dimensionality of MHD flows can be seen from the curl of the Navier-Stokes equation and the cur l of Ohm's law within the quasi-static approximation, which read:

∂ z J z = -∇ ⊥ • J ⊥ = ρ B (u • ∇ω z + ω • ∇u z ) - ρν B ∇ 2 ω z , (1)
∇ × J = σB∂ z u. (2)
These equations express that horizontal layers of fluid can only be fed by electric current if either inertia or viscous friction exist there to balance the Lorentz force. Furthermore, any electric current "leak" pulled into the core by either of these mechanisms results in the presence of velocity gradients along the magnetic field, and hence in three-dimensionality. Conversely, quasi-two-dimensionality is only possible if none of these forces acts, and ∂ z u = 0 in the core. Velocity gradients would then still exist in the

I c ≃ 0 I b I t = I b I I U b U t = U b
l z h (c) l z >> h.
Hartmann layers because of viscous friction, and accordingly, the entire horizontally divergent electric current flows there. The amount of current available to flow into horizontal planes is determined by the intensity of the forcing driving turbulence (or equivalently, by the total current directly injected into the flow, if it is electrically driven, as on figure 1). The height l z of the region where the the forcing current is consumed by this mechanism determines the distance over which the Lorentz force is able to diffuse momentum. Depending on whether in a structure of size l ⊥ and typical velocity U (l ⊥ ), the mechanism is predominantly viscous or inertial, l z respectively scales as:

l ν z ∼ l 2 ⊥ h Ha or l (N ) z ∼ l ⊥ N 1/2 , (3)
where Ha = Bh(ρ/σν) 1/2 and N = σB 2 l ⊥ /(ρU ) are the Hartmann number and interaction parameter. The scaling for l (N ) z was first proposed by [START_REF] Sommeria | Why, how and when mhd turbulence becomes two-dimensional[END_REF] and experimentally verified by [START_REF] Pothérat | Why, how and when mhd turbulence at low rm become three-dimensional[END_REF]. The ratio of the diffusion length l z to the height of the channel determines the dimensionality of the flow, leading to three different cases illustrated on figure 1. If l z ≪ h, momentum diffusion by the Lorentz force does not reach out to both boundaries of the channel, electric eddy currents spread over l z ≪ h and the flow is three-dimensional, with regions near one or both channel walls (depending on the geometry of the forcing) where only a weak, residual flow exists. If l z ∼ h, structures extend across the whole channel but three-dimensionality persists in the bulk. If l z ≫ h, the flow is quasi-two dimensional and electric current flows almost exclusively in the Hartmann boundary layers. It is important to notice that l (N ) z is scale-dependent and so strictly speaking, for a turbulent flow to be quasi-two dimensional, all scales have to satisfy l z (l ⊥ ) ≫ h. Hence, an intermediate state exists where large scales are quasi-two-dimensional while smaller scales are three-dimensional [START_REF] Sommeria | Why, how and when mhd turbulence becomes two-dimensional[END_REF][START_REF] Klein | Appearance of three-dimensionality in wall bounded mhd flows[END_REF]. The three different regimes of flow dimensionality correspond to different electric current paths and therefore different levels of Joule dissipation. We shall now characterise them through the relation between forcing and the measured Reynolds number. For this, we start by noticing that the total current induced by the forcing, or directly injected into the flow, I spreads between the bulk (I c ), and through each of the Hartmann layer (I b and I t ), so that I = I c + I b + I t . Since I c diverges into the core over the height of the structure,

I c ∼ 2πl ⊥ max{l z , h}J c
⊥ . The horizontal current density J c ⊥ is then estimated from (1) using only either the first or the second term in the RHS of the equation, depending on whether viscous effects or inertia dominate, respectively. Using Re b = U b l ⊥ /ν and Re 0 = (I/2πν √ ρσ) to express these scalings non dimensionally (indices b and t refer to bottom and top plates), it comes that if l z ≫ 1 then the flow is quasi-two-dimensional [START_REF] Sommeria | Experimental study of the two-dimensionnal inverse energy cascade in a square box[END_REF][START_REF] Alboussière | Quasi-2d mhd turbulent shear layers[END_REF] and 4), ( 5) and ( 6).

Re b ∼ Re 0 . ( 4
)
If l z ≪ 1 or l z ∼ 1, the core current is pulled in by inertia [START_REF] Pothérat | Why, how and when mhd turbulence at low rm become three-dimensional[END_REF], and

Re b ∼ Re 3/2 0 . (5)
Furthermore, (4) remains valid for any value of l z /h if the current in the core is pulled by viscous effects.

Regimes where this scaling holds shall therefore be called inertialess.

While scalings near the bottom wall give a measure of turbulent intensity where it is forced, scalings near the top wall give a measure of its intensity away from where the forcing is applied. It turns out that if l z ∼ h or l z ≫ h, then scalings for Re t are essentially the same as for Re b , albeit for a small correcting factor [START_REF] Pothérat | Why, how and when mhd turbulence at low rm become three-dimensional[END_REF]. If l z ≪ h, on the other hand, the top wall is outside of the region where turbulence diffuses under the action of the Lorentz force. Any residual flow there is viscously entrained by the neighbouring turbulence, and damped by the Lorentz force. This balance provides a scaling of the form

Re t ∼ Re 1/2 0 , (6)
which characterises the residual flow in this region. All three scalings are observed experimentally to a great precision, over a wide range of values of Re 0 and Ha (see figure 2). For the purpose of the experiment, Re b and Re t were built on the RMS of velocity fluctuations and half of the scale at which energy was injected into the flow ( materialised by the spacing between current injection electrodes L i ).

For low forcing, the inertialess regime dominates both near the wall where turbulence is forced and the wall where it isn't. For high magnetic fields (hence l (N ) z (l ⊥ )/h ≫ 1), a transition takes places to the inertial regime where (5) holds in both regions. At lower fields, where l (N ) z (l ⊥ )/h < 1, this transition is only visible where turbulence is forced. The large scales of turbulence do not reach the top wall and (6) becomes valid in this region.

Can magnetic fields enhance turbulence ?

An important feature of the measurements of turbulent intensity on figure 2 is that the prefactor in all scalings laws linking Re b to Re 0 depends on Ha. This wasn't foreseen in theory and must be attributed to the fact that while the injected current I induces an average flow, we diagnosed turbulent fluctuations, which only receive a fraction of the energy of the average flow. It appears that the ratio of the RMS of velocity fluctuations to the average velocity not only depends on Ha but monotonically increases with it. In other words, the relative turbulent intensity increases with the magnetic field, in apparent contradiction to the widely accepted idea that the magnetic field should suppress turbulence ! The paradox can be resolved by considering the energy balance integrated across the channel and observing that most of the power injected by the forcing P is dissipated ohmically. Using the expression of the Lorentz force put forward by [START_REF] Sommeria | Why, how and when mhd turbulence becomes two-dimensional[END_REF], ohmic dissipation integrated over the channel is expressed as:

h 0 ǫ J dz = -σB 2 h 0 u • ∆ -1 ∂ zz u dz ∼ -σB 2 l 2 ⊥ h u • ∂ 2 zz u + u ′ • ∂ 2 zz u ′ ,
where brackets and prime respectively refer to average and fluctuating quantities. Approximating the derivatives as

∂ 2 zz u ∼ U b 2 /l 2 z and ∂ 2 zz u ′ ∼ U ′ b /l ′2
z , and using (3) with reference velocities U b and U ′2 b 1/2 to evaluate respectively l z and l ′ z , it comes that:

h 0 ǫ J dz ∼ -ρ h l ⊥ U b 3 1 + α 3 , (7)
where α = U ′2 b 1/2 / U b is the relative intensity of turbulent fluctuations. The total power driving the flow across the channel is determined by the injected current as:

h 0 P I dz ∼ IB 2πl ⊥ U b , (8)
and by virtue of (5), equating ( 8) and ( 7) for Ha >> 1, leads to an estimate for α:

α ∼ λHa 1/3 Re -1/6 b . (9)
Experimental data support the scaling Ha 1/3 and confirm that fluctuations become more intense as the field is increased. It must be noted that this scaling applies to fully developed turbulence, for which α 1. The underlying mechanism is that while the magnetic field does suppress turbulence, it does so by elongating structures. In the process, velocity gradients are reduced, and eddy currents responsible for the ohmic dissipation are suppressed. The resulting structures are closer quasi-two-dimensionality and dissipate much less energy than their three-dimensional counterpart. As a result, turbulence retains more kinetic energy. In this sense, the application of a static magnetic field can indeed enhance turbulence.

Figure 1 :

 1 Figure 1: Schematic representation of generic flow configurations in a channel in an external magnetic field. For illustration purposes, the flow is driven electrically by injecting electric current at the bottom wall. Hartmann layers are represented in light blue, paths of electric current in red, and fluid flow in blue. (a) l z << h (b) l z h (c) l z >> h.

Figure 2 :

 2 Figure 2: Re ′b Ha -1/3 and Re ′ t Ha -1/3 vs. Re 0 , representing the RMS of velocity fluctuations near the bottom wall. Solid, dashed-dotted and dashed lines respectively represent scaling laws (4), (5) and (6).