Simon Jones

Mathias Legrand

On Solving One-Dimensional Partial Differential Equations With Spatially Dependent Variables Using The Wavelet-Galerkin Method

Keywords:

The discrete orthogonal wavelet-Galerkin method is illustrated as an effective method for solving partial differential equations (PDE's) with spatially varying parameters on a bounded interval. Daubechies scaling functions provide a concise but adaptable set of basis functions and allow for implementation of varied loading and boundary conditions. These basis functions can also effectively describe C 0 continuous parameter spatial dependence on bounded domains. Doing so allows the PDE to be discretized as a set of linear equations composed of known inner products which can be stored for efficient parametric analyses. Solution schemes for both free and forced PDE's are developed; natural frequencies, mode shapes, and frequency response functions for an Euler-Bernoulli beam with piecewise varying thickness are calculated. The wavelet-Galerkin approach is shown to converge to the first four natural frequencies at a rate greater than that of the linear finite element approach; mode shapes and frequency response functions converge similarly.

Introduction

Mathematical modeling of physical systems often requires the use of a governing partial differential equation (PDE). Phenomena such as elastodynamics of a beam, heat conduction through a plate, or fluid flow through a pipe can all be simplified to respective PDE's with known analytic solutions. When the system of interest has greater complexity it is still possible to describe the system dynamics using a PDE, though complex geometries can make analytical solutions intractable. In such cases solutions can be approximated using numerical methods, such as the finite element method, the finite difference method, the boundary element method, method of lines, general Galerkin methods, or multigrid methods [START_REF] Gockenbach | Partial Differential Equations: Analytical and Numerical Methods[END_REF][START_REF] Fletcher | Computational Galerkin Methods[END_REF][START_REF] Lapidus | Numerical Solution of Partial Differential Equations in Science and Engineering[END_REF].

Relatively recently a new method of solving PDE's has emerged, known as the wavelet-Galerkin method [START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Strang | Wavelets and Filter Banks[END_REF]. Wavelets are well localized, oscillatory functions which have been shown to provide a basis on unbounded domains L 2 .R/ and bounded domains L 2 OEa; b [START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF]. Unlike the Fourier basis, a wavelet basis can define a sparse representation of piecewise signals, including transients and singularities [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. It is hypothesized that these localized properties of the wavelets make them useful test and weighting functions for use in bounded domain Galerkin approaches where problem parameters or field variables feature high gradients or even C 0 continuity.

A number of wavelet functions have been developed to accentuate desirable characteristics, such as orthogonality, compact support, symmetry, and number of vanishing moments. References such as Mallat [START_REF] Mallat | A wavelet tour of signal processing[END_REF], Strang and Nguyen [START_REF] Strang | Wavelets and Filter Banks[END_REF], and Daubechies [START_REF] Daubechies | Ten lectures on wavelets[END_REF] provide excellent detail on the construction and theory behind a wide range of popular wavelets.

The plethora of wavelets from which to select has resulted in extensive literature regarding the wavelet-Galerkin method. For instance, Williams and Amaratunga provide a review of orthogonal wavelet use in engineering [START_REF] Williams | Introduction to wavelets in engineering[END_REF] and specifically to solutions of linear boundary value problems [START_REF] Amaratunga | Wavelet-Galerkin solution of boundary value problems[END_REF]. Ala et al. [START_REF] Ala | An advanced numerical model in solving thin-wire integral equations by using semi-orthogonal compactly supported spline wavelets[END_REF] use semiorthogonal spline wavelets to solve the linear electric field integral, while Joly et al. [START_REF] Joly | Towards a method for solving partical differential equations by using wavelet packet bases[END_REF] use wavelet packets to define adaptive methods of solving linear time dependent PDE's. The wavelet-Galerkin method has also been used for nonlinear PDE's, for example by Beylkin and Keiser [START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF] to efficiently capture shocklike responses in equations described by the semigroup approach. Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] provide a seminal derivation of the discrete orthogonal wavelet-Galerkin method on a bounded interval while investigating the solution to Burger's equation, and Lakestani et al. [START_REF] Lakestani | Numerical solution of the controlled Duffing oscillator by semi-orthogonal spline wavelets[END_REF][START_REF] Lakestani | The solution of a second-order nonlinear differential equation with Neumann boundary conditions using semi-orthogonal B-spline wavelets[END_REF] consider Duffing's equation with various boundary conditions using B-spline wavelets. Restrepo and Leaf [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF] look specifically at periodic solutions using orthogonal wavelets. Periodic wavelets are not generally scale-invariant due to "wraparound" thus special consideration is required when changing the wavelet-scale during the investigation (scale-invariance requires wavelets at any scale to be a pure dilation of the mother-wavelet). Pernot and Lamarque [START_REF] Pernot | A wavelet-Galerkin procedure to investigate time-periodic systems: transient vibration and stability analysis[END_REF] also investigate time-periodic systems, focusing on transient vibrations and stability analysis. Much work has also been done regarding wavelet-Galerkin solutions on bounded domains; the standard wavelet-Galerkin formulation requires the use of wavelets which span outside the domain boundaries thus amendments must be made to the functions to prevent boundary issues. Beylkin [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF], Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], and Romine and Peyton [START_REF] Romine | Computing connection coefficients of compactly supported wavelets on bounded intervals[END_REF] all investigate the computation of inner products and other operators of discrete orthogonal wavelets on bounded domains, while Plonka et al. [START_REF] Plonka | On the construction of wavelets on a bounded interval[END_REF] account for domain bounds using the Chebyshev transform.

In all the work referenced above, the differential equations being solved involved constant parameters. Recently Hein and Feklistova [START_REF] Hein | Free vibrations of non-uniform and axially functionally graded beams using Haar wavelets[END_REF] consider a wavelet-Galerkin approach to solving free-vibration problems involving beams with smoothly varying properties using Haar wavelets. Solutions to bounded PDE's with spatially dependent variables can be difficult to approximate accurately using standard Fourier-Galerkin methods; if parameters differ at the domain boundaries this will act as a jump discontinuity for the periodic Fourier functions and slow solution convergence. Other non-periodic polynomial basis (i.e. Legendre, Chebyshev) can be used effectively for such problems [START_REF] Fletcher | Computational Galerkin Methods[END_REF][START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF], but if parameter variation is only C 0 continuous, a large number of polynomials will be required to capture the rapid gradient changes which again slows convergence. The problem could also be spatially discretized using a finite element approach, but the mesh density must be carefully selected to accurately capture any property variation. In such cases, it is hypothesized that wavelets may allow for a more efficient representation of parameter variation.

In the current investigation, a method for solving one-dimensional partial differential equations with C 0 continuous, spatially dependent variables is illustrated using Daubechies wavelets [START_REF] Daubechies | Ten lectures on wavelets[END_REF]. Derivations for various boundary and loading conditions are given, and solution schemes for both homogeneous and forced equations are detailed.

The paper is broken into three sections. In Section 2 the notation and relevant theory for orthogonal wavelets is reviewed. In Section 3 the discrete orthogonal wavelet-Galerkin method is de-tailed using transverse vibrations of a Euler-Bernoulli beam as an illustrative example. Finally, in Section 4 the method is employed for a clamped-pinned beam with a C 0 continuous, spatially varying cross-section; the results for both free and forced vibration are compared to a finite-element approximation to assess the effectiveness of the method.

Notation of Discrete Orthogonal Wavelets

The family of compactly supported orthogonal wavelets were first reported by Daubechies [START_REF] Daubechies | Ten lectures on wavelets[END_REF], and have since been utilized in a number of areas including signal analysis, image compression, and numerical analysis [START_REF] Strang | Wavelets and Filter Banks[END_REF]. As notation for orthogonal wavelets varies in the literature, a brief review is included below for clarity.

The discrete orthogonal wavelet family is defined by a set of L filter coefficients fp `W `D 0; 1; : : : ; L 1g, where L is an even integer. The fundamental two-scale equations in wavelet theory are defined as

.x/ D L 1 X `D0 p ` .2x `/ (1)
.

x/ D L 1 X `D0 . 1/ `p1 ` .2x `/ (2)
where .x/ is the scaling function and .x/ is the wavelet function, with fundamental support over the finite intervals OE0; L 1 and OE1 L=2; L=2, respectively. These equations can be used to determine the value of the scaling and wavelet function at dyadic points x D n=2 J , n D 0; 1; : : : using the algorithm provided by Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF].

The wavelet filter coefficients p `were derived by Daubechies to produce scaling and wavelet functions with specific properties [START_REF] Daubechies | Ten lectures on wavelets[END_REF][START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], some of which include: the area under the scaling function is unity

Z 1 1 .x/ dx D 1 (3)
the coefficients sum to two

L 1 X `D0 p `D 2 (4)
the scaling function and its translates are orthogonal

Z 1 1 .x/ .x k/ dx D ı 0;k ; k 2 Z (5)
the scaling and wavelet functions are orthogonal

Z 1 1 .x/ .x k/ dx D 0; k 2 Z (6)
the wavelet function has L=2 vanishing moments

Z 1 1
x m .x/ dx D 0; m D 0; 1; : : : ; L=2:

A basis of L 2 .R/ can be formed from the linear spans of the scaling and wavelet functions at resolution level J , respectively [START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF]

J;k .x/ D 2 J =2 .2 J x k/ k 2 Z (8) J;k .x/ D 2 J =2 .2 J x k/ k 2 Z: (9)
The span of scaling functions at level J is commonly denoted V J , while the wavelet span is denoted W j . The characteristics of the wavelet basis allows for a multiresolution analysis of a bounded domain on L 2 OEa; b by decomposition of the space into a chain of closed subspaces with care taken to properly account for functions at the boundaries [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF] V

0 V 1 V 2 V 3 L 2 OEa; b (10)
such that lim

J !1 V J D L 2 OEa; b: (11)
Noting that W J is the orthogonal complement to

V J in V J C1 V J C1 D V J ˚WJ (12)
Equation ( 10) implies

V J C1 D V 0 ˚W0 ˚W1 ˚ ˚WJ : (13)

Discrete Orthogonal Wavelet-Galerkin Method

One arguable drawback to orthogonal wavelets is the two-scale functional form of the fundamental equations ( 1) and ( 2). The corresponding scaling and wavelet functions are highly non-smooth and fractal in nature: as one increases the resolution of either function the shape does not converge but rather continues to increase in complexity [START_REF] Welstead | Fractal and wavelet image compression techniques[END_REF]. This makes accurately estimating inner products of scaling or wavelet functions using numerical integration or quadrature error prone. Interestingly, the exact solution to many orthogonal wavelet operators have been derived using the two-scale relation allowing accurate implementation of the wavelet-Galerkin scheme [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF][START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Romine | Computing connection coefficients of compactly supported wavelets on bounded intervals[END_REF][START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF]. The equation of motion for a Euler-Bernoulli beam is used below as an illustrative example of such an implementation.

Transverse Vibrations of a Euler-Bernoulli Beam

The equation of motion for transverse, damped vibrations of a clamped-pinned Euler-Bernoulli beam assuming no body forces is [START_REF] Doyle | Wave Propagation in Structures[END_REF] E.x/I.

Here y is the transverse beam displacement, E is the elastic modulus, I is the second moment of area, is the mass per unit volume, A is the cross-sectional area, L b is the length of the beam, Á is the damping factor, ! is the frequency of vibration and i D p 1. The beam in this investigation is assumed rectangular in crosssection with a base dimension b and height h. As shown in Figure 1, the variation in height with respect to the x-direction is only C 0 continuous, due to its piece-wise construction. This implies both I D 1=12bh 3 and A D bh are also C 0 continuous functions of x; all other parameters are assumed to remain constant. Accounting for the assumed harmonic response and the spatial dependency of u.x/, I.x/ and A.x/ results in the following equation of motion

ˇ.!/ D ! 2 i !Á E : (17)
The total length of the beam is L b D 2m. Let the domain x 2 OE0; L b be discretized by 2 J points, hence the distance between discrete points x is given by x D L b =2 J . Furthermore, it is convenient in the following sections to consider a domain of OE0; 1, thus a change of variable is introduced

x D x=L b : (18)

Galerkin Approximation

In accordance with the Galerkin method [START_REF] Fletcher | Computational Galerkin Methods[END_REF], an approximate trial solution for u.x/ is introduced using scaling functions at level J as the test functions

u.x/ Š 2 J 1 X kD2 L u k J;k .x/ D 2 J =2 2 J 1 X kD2 L u k .2 J x k/ (19)
where J > 0. The limits of summation are selected such that an approximation of u.x/ can be made for the whole domain x 2 OE0; 1. When x D 0, the summation in Equation ( 19) involves scaling functions evaluated at integer values k. The scaling function .z/ has only L 2 non-zero integer values corresponding to z D 1; 2; : : : ; L 2 [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]. Therefore, the only non-zero J;k .x/ correspond to k D 2 L; 3 L; : : : ; 1.

Similarly when x D 1 the only non-zero values correspond to k D 2 J L C 2; 2 J L C 3; : : : ; 2 J 1. Therefore the required summation spans the values k D 2 L; 3 L; : : : ; 2 J 1.

It is important to stress that Equation ( 19) is a reduced basis approximation of the function u.x/ because no wavelets are included. According to Equation (13) the exact function can be represented as

u.x/ D X k u k J;k .x/ C N X j DJ X `wj;` j;`. x/ ( 20
)
where N is the maximum resolution given by the sampling rate of the function u.x/ [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. As discrete orthogonal wavelet functions are analogous to high-pass filters [START_REF] Strang | Wavelets and Filter Banks[END_REF], neglecting the wavelet terms in Equation ( 20) is akin to removing portions of the high frequency content of u.x/, resulting in a smoothed approximation of the true function. As shown in Equation (10) however, the scaling function approximation approaches L 2 OEa; b as the level J is increased. Therefore to improve the accuracy of the approximation the wavelet functions at level J can be included, or more simply the level of the scaling functions can be increased to J C 1 as detailed in Equation [START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF]. Restricting the test functions to scaling functions at a single level simplifies the algorithms while also providing a reduced basis which is adaptable to the desired level of accuracy.

Another consideration before continuing with the Galerkin approach is how to deal with the spatially dependent terms I.x/ and A.x/ from Equation ( 16). Since the height function h.x/ is piecewise linear, directly substituting this dependence into the governing equation would result in different equations to be satisfied depending on the x location. If one considers a more general case where the height dependence is not piecewise linear, the Galerkin formulation would require finding the integral of a triple product of two scaling functions and an arbitrary function of x (e.g. the integral of the weighting function, test function, and the I.x/ function). There are a limited set of functions where this triple product integral can be found exactly [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]; in general this integral would have to be approximated using numerical integration which is known to be error prone [START_REF] Welstead | Fractal and wavelet image compression techniques[END_REF].

The novelty of the proposed approach lies in the decomposition of the spatially varying functions, A.x/ and I.x/, using discrete orthogonal scaling functions such that

A.x/ Š 2 J 1 X j D2 L A j J;j .x/ D 2 J =2 2 J 1 X j D2 L A j .2 J x j / (21) I.x/ Š 2 J 1 X j D2 L I j J;j .x/ D 2 J =2 2 J 1 X j D2 L I j .2 J x j /: (22)
The scaling function approximations provide relatively accurate representations of the piecewise functions due to the vanishing moment properties in Equation ( 7) allowing perfect reproduction of polynomials up to order .L=2 1/ [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]. This approach results in a single governing equation, regardless of x location, and the resulting triple product integrals can be found exactly (see below for full explanation).

The coefficients A j and I j are calculated using the inner product [START_REF] Mallat | A wavelet tour of signal processing[END_REF] A j D Z 1 0 A.x/ J;j .x/ dx and I j D Z 1 0 I.x/ J;j .x/ dx: [START_REF] Welstead | Fractal and wavelet image compression techniques[END_REF] This inner product can be performed using the discrete wavelet transform (DWT) or fast wavelet transform (FWT) [START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Strang | Wavelets and Filter Banks[END_REF]. Attention must be paid to the calculation of the coefficients near the boundary x D 0. Standard decomposition requires the scaling functions to extend into the negative x domain by L 2 points, thus a decision on how to represent this boundary must be made. A number of options exist including scaling function extrapolation [START_REF] Amaratunga | Wavelet-Galerkin solution of boundary value problems[END_REF], symmetric extension [START_REF] Mallat | A wavelet tour of signal processing[END_REF], or the use of boundary scaling functions [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF][START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Romine | Computing connection coefficients of compactly supported wavelets on bounded intervals[END_REF]]. In the current investigation it is assumed a symmetric extension is appropriate (i.e. A. x/ D A.x/); due to the beam geometry in this problem all the methods listed would actually provide the same results.

The Galerkin method defines the residual R as the substitution of Equations ( 19), ( 21) and ( 22) into Equation ( 16)

R D 2 J 1 X j D2 L 2 J 1 X kD2 L h I j u k .0/ J;j .x/ .4/ J;k .x/ C 2 .1/ J;j .x/ .3/ J;k .x/ C .2/ J;j .x/ .2/ J;k .x/ ˇ.!/A j u k .0/ J;j .x/ .0/ J;k .x/ i ( 24
)
where the superscript (n) refers to the order of differentiation.

The nth derivative of the scaling function .n/ J;k .x/ is defined as

.0/ J;k .x/ D J;k .x/ (25)
.n/

J;k .x/ D d n J;k dx n .x/ D d dx .n 1/ J;k .x/ (26)
By amalgamating Equations ( 8), ( 18) and ( 26), and accounting for the vanishing moment condition in Equation [START_REF] Daubechies | Ten lectures on wavelets[END_REF], it is possible to write

.n/ J;k .x/ D 2 J =2
x n .n/ .2 J x k/; n D 0; 1; : : : ; L=2 1 [START_REF] Humar | Dynamics of structures[END_REF] It is important to note the limit of n D L=2 1 due to the vanishing moment condition; in the current investigation where a fourthorder derivative is required this condition implies a minimum scaling function order of L D 10. A method for determining values of the scaling function and its derivatives at dyadic points x D k=2 J is given by Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF].

Continuing with the Galerkin method, scaling functions of level J are selected as the weighting functions and the inner product of the residual and the weighting functions is set to zero [START_REF] Fletcher | Computational Galerkin Methods[END_REF] Z 1 0 R J;`. x/ dx D 0; `D 2 L; 3 L; : : : ; 2 J 1 (28) which, using Equation [START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF], results in

2 J 1 X j D2 L 2 J 1 X kD2 L I j u k a `;j;k C 2b `;j;k C c `;j;k ˇ.!/A j u k .
Chen et al. refer to triple product integrals of this form as threeterm connection coefficients [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] and show that the exact solution to these proper integrals can be found. Reference [START_REF] Zhang | Comments on 'The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] suggests improved algorithms for calculation of these connection coefficients. The standard notation for the three-term connection coefficients is

H `;k D 2 J =2 2 J 1 X j D2 L A j 0;0 j `;k `.2 J `/ 0;0 j `;k `.
`/Á for k; `D 2 L; 3 L; : : : ; 2 J 1:

It should be noted that the three-term connection coefficients are parameter independent, thus they can be stored in matrix form and reused efficiently for parametric analyses or even different PDE's requiring the same scaling function derivative combinations. This makes the wavelet-Galerkin method an efficient tool overall.

Imposing Boundary Conditions

The final step before a solution can be found involves prescription of the appropriate boundary conditions on the basis functions. As shown in Figure 1, the beam in this investigation is clamped at x D 0 and pinned at x D 1, which is equivalent to u.0/ D 0, u ;x .0/ D 0, u.1/ D 0, and E O I u ;xx .1/ D 0 respectively, where O I is the moment of inertia at x D 1. Consider the clamped boundary conditions at x D 0; using Equation ( 19) these can be written in the scaling function domain as

u.0/ D 2 J =2 2 J 1 X kD2 L u k . k/ D 0 (41) u ;x .0/ D 2 J =2 x 2 J 1 X kD2 L u k .1/ . k/ D 0 (42)
Noting that .x/ and .1/ .x/ have non-zero values for integers x D 1; 2; : : : ; L 2 only [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], these equations simplify to

1 X kD2 L u k . k/ D 0 (43) 1 X kD2 L u k .1/ . k/ D 0: (44)
Now consider the pinned boundary condition at x D 1. If the end is free of any applied moments, the BC's can be written

u.1/ D 2 J =2 2 J 1 X kD2 L u k .2 J k/ D 0 (45) E O I u ;xx .1/ D 2 J =2 x 2 E O I 2 J 1 X kD2 L u k .2/ .2 J k/ D 0: (46)
Accounting again for the non-zero values of .x/ and .2/ .x/, these can be simplified to

2 J 1 X kD2 J LC2 u k .2 J k/ D 0 (47) 2 J 1 X kD2 J LC2 u k .2/ .2 J k/ D 0: (48)
If instead there is a moment M.t / D M e i !t applied at x D 1, as depicted in Figure 1, the respective BC becomes

2 J 1 X kD2 J LC2 u k .2/ .2 J k/ D Mx 2 2 J =2 E O I : (49)

Free Vibrations Solution Scheme

The natural frequencies of the beam can be found by assessing the homogeneous solution to Equation (39). This requires setting the forcing moment M D 0 and the damping factor Á D 0. Equation ( 39) is a generalized eigenvalue problem in ˇ.!/: the eigenvalues can be used to determine the natural frequencies of the system ! n using Equation ( 17); the corresponding eigenvectors can be used to compute the modeshapes using Equation [START_REF] Romine | Computing connection coefficients of compactly supported wavelets on bounded intervals[END_REF]. Thus the four boundary conditions (43), (44), ( 47) and (48) must be imposed in a manner so as not to disturb the generalized eigenvalue form.

One efficient approach is to constrain four participation factors using the boundary conditions. To accomplish this it is convenient to first reorder Equation (39) so that

O G ˇ.!/ O H u m k u s k D 0 (50)
where u m k are the master scaling function participation factors and u s k are the slave factors which will be constrained. O G and O H are the reordered G and H matrices, respectively. The size of u m k is .2 J C L 6/ 1 and u s k is 4 1 to account for the four boundary conditions.

Boundary conditions (43), ( 44), ( 47) and ( 48) can be used to form the constraint equations

u m k u s k D O Cu m k (51)
where O C > D OEI > ; C > is the constraint matrix. I is an identity matrix of order .2 J C L 6/ and C is a 4 .2 J C L 6/ matrix derived from Equations ( 43), ( 44), ( 47) and (48). The slave participation factors should be selected such that the condition number of C is minimized; for the current boundary conditions the slave factors are u 1 , u 2 , u 2 J 2 , u 2 J 1 . A sample calculation of C is given in the Appendix.

Substituting Equation (51) into (50) and pre-multiplying by the transpose of O C results in the constrained system of equations

O C > O G ˇ.!/ O H O Cu m k D 0 (52)
thus the boundary conditions have been imposed without disturbing the generalized eigenvalue form of the equation. The eigenvector coefficients for the slave participation factors can be determined using Equation (51).

Forced Vibrations Solution Scheme

As an example of how loading can be incorporated in the model, a harmonic moment loading condition at the right end of the beam is considered as shown in Figure 1. It should be noted that more general loading conditions, such as distributed loads, can also be included using the wavelet-Galerkin method. Similar to the free response case, the system of equations (39) is constrained using the appropriate boundary condition (43), (44), (47) and (49). Due to the applied moment in Equation (49) the constraint equation must be amended to

u m k u s k D O Cu m k C O F: (53)
O F is of size .2 J C L 2/ 1 and all entries are zero except for the row associated with Equation (49) in C where the entry is .Mx 2 /=.2 J =2 E O I /. The resulting constrained system of equations is

O C > O G ˇ.!/ O H O Cu m k D O C > O G ˇ.!/ O H O F ( 54
)
This system can be solved using standard matrix inversion for participation factors u m k . The solution y.x; t / for the beam can be found by solving for u.x/ using Equation ( 19), finding y.x; t / D u.x/e i !t , and finally transforming back to the x domain using x D L b x.

Analysis and Discussion

The discrete orthogonal wavelet-Galerkin method is employed to approximate the transverse vibrations of an Euler-Bernoulli beam as described by Equation [START_REF] Lakestani | Numerical solution of the controlled Duffing oscillator by semi-orthogonal spline wavelets[END_REF]. The beam, whose properties are given in Table 1, varies in thickness as shown in Figure 1. The 10 coefficient Daubechies scaling function (DB10) is used for all investigations; the level J is adjusted to determine the effect on prediction accuracy. A finite element model (FEM) comprised of linear beam elements [START_REF] Humar | Dynamics of structures[END_REF] is also constructed using the beam parameters listed in Table 1 to act as a comparison model.

Free Vibration

Figure 2 shows the convergence of the first four natural frequencies using the finite element model (FEM) and the wavelet-Galerkin model (W-G); the number of unknowns shown equate to J D 5; 6; 7; 8. As shown in the figure, the wavelet-Galerkin approach converges to a solution more rapidly than the finite element method. This is attributed to the scaling functions' ability to better represent the piecewise variation in cross-sectional area and second moment of inertia than can the linear finite elements. The FEM discretization results in a step-wise approximation of the thickness variation thus a fine mesh is required to converge to a solution. The vanishing moment properties of the Daubechies scaling functions allow perfect reconstruction of the linear variations in area and cubic variation in moment of inertia with only localized error at the piecewise joints. These joint locations have discontinuous slope and require higher level wavelets to perfectly reconstruct, hence the convergence as J is increased.

For this example, computational speed is essentially proportional to the number of unknowns. The element matrices can be pre-computed and stored for the FEA method in a similar way to the connection coefficients for the W-G method. The matrices of interest are relatively small thus the time required for the preand post-multiplication in Equation ( 52) and ( 54) is negligible. The computational time required for the FEA or W-G method is thus governed by the solution to the eigenvalue or linear algebra problem, which is the same for both methods. It is currently unclear how the wavelet-Galerkin method would scale up to much larger problems; this is an active area of research.

The associated normalized displacement mode shapes predicted by FEM and W-G at J D 7 are shown in Figure 3. The predicted mode shapes from both models overlay each other (< 0:1 % relative error rms), attesting to the validity of the purposed wavelet-Galerkin method.

Forced Vibration

The frequency response function for the displacement of the midspan of the beam (x D 1 m) due to a 10 N m harmonic moment at x D 2 m is provided in Figure 4. Included in the figure is the FEM prediction for 128 degrees of freedom and the W-G prediction at levels J D 4; 5; 6. The peaks of the curve correspond to the natural frequencies shown in Figure 2; these peaks are bounded due to the damping included for this portion of the investigation (see Table 1).

The plot shows how the wavelet-Galerkin solution converges rapidly as the level J is increased. As in the free vibration case, the values of the natural frequencies are over predicted at smaller J , most noticeably for higher modes. This is attributed to an inadequate number of scaling functions necessary to precisely reconstruct the increasingly rapid changes in slope of the higher modes, as well as the thickness variation of the beam.

Conclusion

The discrete orthogonal wavelet-Galerkin method is introduced as an efficient method of solving partial differential equation (PDE) with spatially varying parameters on a bounded domain. The Daubechies scaling functions provide a concise but adaptable set of basis functions: increasing the level of the scaling functions results in a more encompassing subspace of L 2 OEa; b and thus a more accurate approximation to the solution of the PDE. The spatial dependence of the parameters can be described using the same scaling functions, which allows the PDE to be written as a set of linear equations composed of three-term connection coefficients. The connection coefficient matrices are parameter and even problem independent thus they can be stored and reused for efficient parametric analyses of various PDE's.

An example problem of transverse vibrations of a pinnedclamped Euler-Bernoulli beam with piecewise varying thickness is used to describe the method. Solution schemes for both homogeneous and forced PDE's are developed for a number of boundary conditions. The natural frequencies computed using the wavelet-Galerkin method are shown to converge more rapidly than do those predicted by a linear finite element model with equivalent degrees of freedom. The associated modeshape are also shown to match between the two models. Finally, the frequency response function predicted by the wavelet-Galerkin method is shown to converge as the level of scaling function is increased as expected. These results suggest the wavelet-Galerkin method is well suited to approximating the solution of PDE's with spatially dependent parameters. where a 11 D . .0/ .1/ .

Figure 2 .

 2 Figure 2. Convergence of the first four natural frequencies using the finite element method (FEM) and wavelet-Galerkin method (W-G)

Figure 3 .

 3 Figure 3. First four mode shapes predicted by the finite element method (FEM) and wavelet-Galerkin method (W-G) using 128 degrees of freedom (J D 7)

Figure 4 .

 4 Figure 4. Frequency response function at midspan (x D 1 m) under a 10 N m harmonic moment acting at x D 2 m

  Thickness variation h.x/ of the beam where

	y	x	1.00cm	1.50cm	M(t)
	0		0.50 0.75	1.25 1.50	2.00m
			Figure 1.

I.x/u ;xxxx .x/ C 2I ;x .x/u ;xxx .x/ C I ;xx .x/u ;xx .x/ ˇ.!/A.x/u.x/ D 0 (16)

Table 1 .

 1 Beam properties

  1/ .4/ .0/ .4/ .1/ .1//= a (58) a 12 D . .0/ .1/ .1/ .3/ .0/ .3/ .1/ .1//= a (59) a 21 D . .0/ .4/ .1/ .2/ .0/ .2/ .1/ .4//= a (60) a 22 D . .0/ .3/ .1/ .2/ .0/ .2/ .1/ .3//= a (61) b 11 D . .0/ .1/ .2/ .4/ .0/ .4/ .2/ .1//= b (62) b 12 D . .0/ .1/ .2/ .3/ .0/ .3/ .2/ .1//= b (63) b 21 D . .0/ .4/ .2/ .2/ .0/ .2/ .2/ .4//= b (64) b 22 D . .0/ .3/ .2/ .2/ .0/ .2/ .2/ .3//= b (65) and a D .0/ .2/ .1/ .1/ .0/ .1/ .1/ .2/ (66) b D .0/ .2/ .2/ .1/ .0/ .1/ .2/ .2/:

								(67)
	Combining these gives
	8 < :u 2 u 1 u 14 u 15	9 > > = > > ;	D	2 6 6 4 a 11 a 12 0 a 21 a 22 0 0 0 b 11 b 12 0 0 0 0 b 21 b 22	3 7 7 5	8 < :u 4 u 3 u 0 u 1 u 11 u 12 u 13 : : :	9 > > > > > > ; > > > > > > > > > > > > > > > > =	(68)
	or u s k D Cu m k .

Acknowledgements

This work was supported in part by the Tomlinson Postdoctoral Fellowship.

Appendix

Here a sample constraint matrix C from Section 3.4 is constructed using the boundary conditions given in Equations ( 43), (44), (47) and (48) for L D 6 and J D 4. Substituting these parameters into the respective boundary condition equations results in


W J Liu

D Krasnov

J Schumacher

C Karcher

A Thess

DNS of natural convection in liquid metal with strong magnetic fields in rectangular geometry

Keywords:

Direct numerical simulations (DNS) of natural convection in liquid metal within closed boxes heated from below with uniform strong vertical magnetic fields are conducted. The main aim is to explore the possibilities and mechanisms of convection instabilities in such flows, which are referred as MHD convection. It is shown that vertical magnetic field can suppress convection at the center and reside the convection rolls attach to side walls. Stronger magnetic field can even generate new convection rolls which are aligned with the direction of magnetic field.

Introduction

We discuss the interaction between a buoyancy-induced flow and a strong magnetic field. It is known that turbulence can be fully suppressed by strong magnetic fields, leading to a laminar regime. However, this does not mean that laminar steady-state regime is always correct. The results in experiment [START_REF] Genin | 8th Pamir Conference on Fundamental and Applied MHD[END_REF] where the flows through a pipe both with magnetic field and temperature gradient are considered show that temperature fluctuations disappear at moderate magnetic fields, but reappear at stronger magnetic fields in form of high-amplitude low-frequency oscillations indicating that convection-related instabilities may occur. This phenomenon is expected to appear in all systems with strong magnetic fields and buoyancy effects, particularly, in the case of liquid metal blankets of a fusion reactor, where Hartmann and Rayleigh numbers are very high. A detailed investigation of thermal convection in liquid metal flows with strong magnetic fields appears necessary. In this work, we consider the influence of strong magnetic field on natural convection in liquid metal within a rectangular box, which is relevant to the situation with weak or zero imposed flow, for example, in the case of the HCLL blanket, where the liquid metal is used only for breeding. We focus on fundamental analysis rather than on specific applications.

Equations

We consider flows in a closed box with imposed vertical magnetic field and vertical temperature gradient created by maintaining upper and bottom walls at constant and different temperatures, as shown in figure 1. The other four walls are thermally insulated and all walls are perfectly electrically insulating. The fluid is considered as an electrically conducting, incompressible, Newtonian fluid, e.g. a liquid metal. The Boussinesq approximation is applied for the temperature-related buoyancy force. And we assume that the magnetic Reynolds number and the magnetic Prandtl number are both much smaller than one. The non-dimensional governing equations are:

0   u (1)   2 1 Lb p t Gr           u u u u F F (2)   2 2 1 t GrPr         u (3)
Where

Pr z z A t z A t z Nu Ra u T T z        .
The problem is solved by direct numerical simulation with FORTRAN code solver based on the projection-type and nearly fully conservative finite difference scheme described in [START_REF] Krasnov | [END_REF]. In the present study, a modified version for nonperiodic boundary conditions is used and the detailed discussion of the method can be found in [3]. In addition, we use semi-implicit scheme for solving temperature transport equation, where convective terms are computed explicitly while diffusive terms are computed implicitly. The code is also adapted for massive-parallel computers, including hybrid parallelization with MPI and Open MP interfaces. , the convection structures are turbulent. Comparing the results of the different aspect ratios, we find that bounded side walls can suppress the development of convection rolls and have obvious influence on the flow structures. In a long rectangular box, the preferred mode is a finite number of large circulation rolls with axes parallel to the short side, e.g. figure 2(b), where 1: 4 :1  .While in a square box, the convection structures can develop in both horizontal directions with large-scale circulations oscillating around corners, e.g. figure 2(c), where 4 : 4 :1  .

Results

Fully

Uniform magnetic fields in the vertical direction are imposed to investigate the influence of magnetic fields on the system and the magnetic fields are applied at the very beginning. The fully developed convection structures under the influence of magnetic fields with Ha=300 are shown in figure 3. We can see that the magnetic fields change the structures of convection and orientation of convection rolls completely. With vertical magnetic field, the original convection rolls tend to reside near the side walls, while the convective motion at the center is suppressed. 4. We can see that with stronger vertical magnetic field, the convection rolls are more regular and more likely to attach to the side walls. Furthermore, new convection rolls appear at the corners which are alined with the direction of the magnetic field. Once the magnetic field is strong enough, convection in the system can be almost fully suppressed, only leaving very weak flow at the corner.

Conclusions

We have conducted direct numerical simulations of natural convection in liquid metal within closed boxes heated uniformly from below with uniform strong vertical magnetic field. The influences of magnetic field and aspect ratio on the convection structures have been analyzed. The magnetic fields can completely change the structures and orientation of convection rolls. The original convection rolls tend to reside near the side walls under vertical magnetic field, while the convective motions at the center tend to be suppressed and new structures alined magnetic fields appear. Bounded walls can influence the development of convection rolls and flow structures. More detailed and quantitative analysis about the influence of the aspect ratio will be conducted.

Fig. 1 :

 1 Fig.1: Sketch of geometry and reference system.

  velocity is no slip boundary condition, that is, the velocity is set zero on all the six walls. For temperature, one pair of the opposite walls (e.g. zz-walls) is set as isothermal walls with constant temperature (e.g.

Fig. 2 :

 2 Fig. 2: Convection patterns without magnetic field for different aspect ratios. Instantaneous temperature field (upper row) and streamtraces (lower row) are showed. Arrows indicate the direction of gravity acceleration and temperature gradient.

Fig. 3 :

 3 Fig. 3: Convection patterns with magnetic field for different aspect ratios. Instantaneous temperature field (upper row) and streamtraces (lower row) are showed. Direction of magnetic field is additionally indicated.

  Fig.4: convection patterns in cubic boxes under vertical magnetic fields with different Ha. Instantaneous temperature field (upper row), 3D streamtraces (middle row) and top view of streamtaces (lower row) are showed.

  The dimensionless control parameters are Prandtl number, Grashof number or Rayleigh number, Hartmann number, and aspect ratio, which are defined by

	). All other walls are adiabatic with 3 2 () hc g T T H Ra Gr Pr     ; 0 Ha B H  In our case there is no imposed flow, the velocity is scaled by freefall velocity 0 n     (e.g. x=±1). Pr    ;   ; ::     y  . Lx Ly H f () 0 at y=±1; hc U   and length is 0     at x (4) g T T H scaled by the height of the geometry H . Time is consequently scaled by the freefall time ff  . Furthermore, the T H U global momentum transfer in the system can be described by Reynolds number rms f  , Re u Re where ff Re U H Ra Pr   is freefall Reynolds number and 2 , rms i V t velocity. The global heat transfer in the system can be described by Nusselt number uu    is the root mean square of the ( ) , ( ) , ( )

  DNS results of natural convection in closed box without magnetic field are shown in figure2. Different aspect ratios are considered. We can see that large-scale circulations exit in the boxes. They are typical of Rayleigh-Bénard convection. And with Ra=10 6

										Ra=10 6 ,
		1:1:1, 1: 4 :1, 4 : 4 :1	and Ha from 100 to 1000 are used in the analysis. And the resolutions are 128×128×128 for
		1:1:1	, 128×512×128 for		1: 4 :1	and 512×512×128 for		4 : 4 :1	.

3D DNS of convection in closed boxes are conducted. The effects of parameter range on the flow structures, i.e. variations in aspect ratio and Hartmann number, are analyzed. Values of Pr=0.025, typical for liquid metals,

Acknowledgment

The authors are grateful to Thomas Boeck and Yurii Kolesnikov for interesting discussions and useful comments. Financial support was provided by the Deutsche Forschungsgemeinschaft (DFG) within the RTG Lorentz Force Velocimetry and Eddy Current Testing. W.J. Liu's work was supported by the China Scholarship Council and D. Krasnov's work was supported by the Helmholtz Alliance "Liquid Metal Technologies". Computer resources were provided by the computing center of TU Ilmenau.