G Bonithon

P Joyot
email: pjoyot@estia.fr

F Chinesta

P Villon

Improving BEM solvers: The Proper Generalized Decomposition Boundary Element Method for solving parabolic problems

Keywords: Boundary element method, Separated representations, Proper Generalized Decomposition

In this work we propose a new approach for solving the heat equation within the Boundary Elements method framework. This technique lies in the use of a separated representation of the unknown field that allows decoupling the space problem (that results steady state) from the temporal one (one dimensional that only involves the time coordinate).

Introduction

The Boundary Elements Method (BEM) allows efficient solution of partial differential equations whose kernel functions are known. The heat equation is one of these candidates when the thermal parameters are assumed constant (linear model). When the model involves large physical domains and time simulation intervals the amount of information that must be stored increases significantly.

We propose here an alternative strategy able to change the nature of the problem. Thus, the temperature field involved by the so called heat equation is approximated using a separated representation involving products of space and time functions. This kind of approximation is not new, in fact proper orthogonal decomposition [START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF] allows such one decomposition, but in this case this decomposition must be applied a posteriori, i.e. on the transient solution of the considered model.

The technique that we propose in this paper allows to transform the transient model in a sequence of space problems (all of them steady state) and time problems (that only involve the time coordinate). This iteration procedure lead to a proper space-time generalized decomposition of the model solution.

The efficiency of such one approach was proven in [START_REF] Ammar | The nanometric and micrometric scales of the structure and mechanics of materials revisited: An introduction to the challenges of fully deterministic numerical descriptions[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part ii: Transient simulation using space-time separated representations[END_REF]. In our knowledge, this technique has never been coupled with a BEM for solving the resulting steady problem defined in the physical domain.

A Proper Generalized Decomposition Boundary Element Method

Let Ω ⊂ R n be a domain of boundary Γ and let τ ⊂ R + a time interval in which the model is defined. The strong form of the heat equation with homogeneous initial and boundary conditions writes: Find the temperature field (assuming known the heat source b (x,t)) such that:

∂ ,t u (x,t) -u (x,t) = b (x,t) x ∈ Ω,t ∈ τ (1) u (x,t) = 0 x ∈ Γ,t ∈ τ (2) u (x, 0) = 0 (3)
In the proposed technique the transient equation is not solved by applying the standard BEM technique. On the contray, we are using the fact that u (x,t) can be expressed into the separated representation form:

u (x,t) ≈ n ∑ k=1 R k (x) S k (t)
We are looking for the functional couples (R k (x), S k (t)). We must recall that the efficiency of the method depends on a number of functional couples n small enough. However, the solution of many models can be written from a reduced number of functional couples.

First, we should give the equations that must verify the different functional couples. For the sake of simplicity we introduce the notation:

R i, j (x) = S i (t) , ∂ ,t S j (t) τ R j (x) -S i (t) , S j (t) τ R j (x) S i, j (t) = R i (x) , R j (x) Ω ∂ ,t S j (t) -R i (x) , R j (x) Ω S j (t)
where

(u, v) τ = τ u(t)v(t)dt and (u, v) Ω = Ω u(x)v(x)dx.
We can prove that Proposition 1 (Galerkin method for separated representations). The couples (R k (x) S k (t)) that verify the strong form of the heat equation are also solution of the non-linear systems:

     R 1,1 (x) + • • • + R 1,n (x) = (S 1 (t) , b (x,t)) τ . . . . . . R n,1 (x) + • • • + R n,n (x) = (S n (t) , b (x,t)) τ      S 1,1 (t) + • • • + S 1,n (t) = (R 1 (x) , b (x,t)) Ω . . . . . . S n,1 (t) + • • • + S n,n (t) = (R n (x) , b (x,t)) Ω (4)
Proof. Let H τ be a Hilbert space with the scalar product

(u(t), v(t)) τ = τ u(t)v(t) dt
We assume that the functions S k (t) define a basis B τ of H τ . By applying the Galerking method to the heat equation, we are looking for the functions

R k (x) such that n ∑ k=1 R k (x) ∂ ,t S k (t) -R k (x) S k (t) , S i (t) τ = 0 ∀i = 1 • • • n
that can be rewritten as:

n ∑ k=1 R i,k (x) = (S i (t) , b (x,t)) τ ∀i = 1 • • • n (5)
Now, we consider another Hilbert space H Ω with the associated scalar product

(u(x), v(x)) Ω = Ω u(x)v(x) dx
We consider that functions R k (x) define the basis B Ω of H Ω . Thus, we look for the functions S k (t)

such that n ∑ k=1 R k (x) ∂ ,t S k (t) -R k (x) S k (t) , R i (x) Ω = 0 ∀i = 1 • • • n
that can be rewritten as:

n ∑ k=1 S i,k (t) = (R i (x) , b (x,t)) Ω ∀i = 1 • • • n (6)
Equations eq [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] and eq (6) correspond to the system eq (4).

In the context of the proposed strategy, the number of functional couples R k (x) S k (t) is "a priori" unknown. We propose to build-up this basis iteratively, until reaching convergence.

Thus, at the first iteration we consider the basis B Ω (resp. B τ ) composed of a single element R 1 (x), (resp. S 1 (t)). Thus, we look for R 1 (x), S 1 (t) such that:

R 1,1 (x) = (S 1 (t) , b (x,t)) τ S 1,1 (t) = (R 1 (x) , b (x,t)) Ω
At iteration p, the basis B Ω (resp B τ ) is composed by the functions previously computed R 1 (x) to R p-1 (x), (resp. S 1 (t) to S p-1 (t)) and by R p (x) and S p (t) verifying

R p,p (x) = S p (t) , b (x,t) τ -∑ p-1 i=1 R p,i (x) S p,p (t) = R p (x) , b (x,t) Ω -∑ p-1 i=1 S p,i (t) (7)
Proposition 2. First, we introduce the notation b

i (x,t) = R i (x) ∂ ,t S i (t) -R i (x) S i (t). If at iteration n, we have b (x,t) = ∑ n k=1 b i (x,t), then u (x,t) = ∑ n k=1 R k (x) S k (t)
is solution of eq [START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF]. We can affirm that the method converges. In that case the couples R k (x)S k (t) verify the system eq (4), as well as the conditions:

∑ n i=p+1 S p (t) , b i (x,t) τ = 0 ∑ n i=p+1 R p (x) , b i (x,t) Ω = 0 1 p n -1
that can be written as:

∑ n i=p+1 R p,i (x) = 0 ∑ n i=p+1 S p,i (t) = 0 1 p n -1 Proof. Because b (x,t) = ∑ n k=1 b i (x,t), we have S p (t) , b (x,t) τ = n ∑ i=1 S p (t) , b i (x,t) τ = n ∑ i=1 R p,i (x)
At iteration p, R p (x) is computed from eq (7) that also writes

R p,p (x) = n ∑ i=1 R p,i (x) - p-1 ∑ i=1 R p,i (x)
This implies ∑ n i=p+1 R p,i (x) = 0.

Remark 1. The proposed iteration scheme stops when the relation

b (x,t) - n ∑ k=1 R i (x) ∂ ,t S i (t) -R i (x) S i (t) Ω×τ ε
holds, where ε is a small enough threshold value.

In order to enforce the homogeneous boundary and initial conditions eq (2), one should enforce

R k (x) = 0 x ∈ Γ and S k (t) = 0 t ∈ τ.
In order to solve the resulting non-linear problem eq (7) at iteration p, we use a fixed point strategy.

We denote by R (q) p (x) and S (q) p (t) the values of R p (x) and S p (t) at iteration q. We consider that R (q) p (x) and S (q) p (t) are known, and we determine R

(q+1) p (x) from S (q) p (t) , ∂ ,t S (q) p (t) τ R (q+1) p (x) -S (q) p (t) , S (q) p (t) τ R (q+1) p (x) = S (q) p (t) , b (x,t) τ -∑ p-1 i=1 R (q) p,i (x) (8)
In fact is that equation related to the boundary conditions on R p (x) that we are solving by applying the BEM. However, as the kernel is different depending on the sign of S (q) p (t) , ∂ ,t S (q) p (t) τ , we rewrite that equation as:

-S (q) p (t) , S (q) p (t) τ R (q+1) p (x) = S (q) p (t) , ∂ ,t S (q) p (t) τ R (q) p (x) + S (q) p (t) , b (x,t) τ -∑ p-1 i=1 R (q) p,i (x) (9)
where the first term of the rigth hand member of eq (8) is computed at iteration q instead of q + 1. Thus, R (q+1) p will be the solution of a Poisson equation.

We are now calculating S

R (q+1) p (x) , R (q+1) p (x) Ω ∂ ,t S (q+1) p (t) -R (q+1) p (x) , R (q+1) p (x) Ω S (q+1) p (t) = R (q+1) p (x) , b (x,t) Ω -∑ p-1 i=1 S (q) p,i (t) (10)
This equation is a simple ODE whose solution does not represent any difficulty.

The method needs the computation of integrals R i (x) , R j (x)

Ω , R i (x) , R j (x) Ω
, and for that we define an approximation in Ω such that R i (x) = N T (x)s, N T (x) is the matrix containing the shape functions and s the vector containing the nodal values. We can chose the finite element shape functions, or the ones associated with the MLS or the RBF. Vector s is obtained using the BEM representation formula.

We proceed in the same way for the integral (S i (t) , S i (t)) τ , (S i (t) , ∂ ,t S i (t)) τ .

Numerical example

In this section we validate the just proposed approach. We consider the square domain Ω = [0, 1] × [0, 1] and the time interval τ = [0, 0.2]. The source term is constant b (x,t) = 1 and both boundary and initial conditions homogeneous.

The exact solution of the problem writes: (see [START_REF] Carslaw | Conduction of Heat in Solids[END_REF])

u as (x,t) = τ 1 0 1 0 G (x, y, x, ȳ, t) d xd ȳd t where G (x, y, x, ȳ, t) = 4 ∞ ∑ n=1 ∞ ∑ m=1 sin (nπx) sin (mπy) sin (nπ x) sin (mπ ȳ) exp -π 2 n 2 + m 2 t
The Poisson's equation eq (9) is solved by the BEM. the domain boundary Γ is discretized by using n Γ × n Γ segments Γ i , the field u (x,t) being constant on each one.

To compute de integrals R i (x) , R j (x)

Ω , R i (x) , R j (x)
Ω the field u (x,t) is reconstructed within Ω on a mesh compatible with the segments defined on the domain boundary Γ i . Each element of the mesh defined a domain Ω i (see figure 1). The field is approximated using quadratic MLS based shape functions.

The integrals are computed by using a single integration point located at the middle point of each segment Ω i . The time interval is discretized using linear MLS based shape functions from the n τ nodes uniformly distributed.

We are checking the convergence for different space (i.e. n Γ ) and time (i.e. n τ ) discretizations.

Table 1 and figure 2 show the evolution of the L 2 error in time and space as a function of the number of functional couples R k (x) S k (t) involved in the approximation of u (x,t) and the mesh densities. We define this error as:

n Γ n τ e 2 1 e 2
e n = n ∑ k=1 R k (x) S k (t) -u as (x,t) L 2
Ω×τ we can notice that for a number of couples n the error e n decreases reaching a plateau when n Γ increases. For reducing this error we must increase n Γ and recompute more couples R k (x) S k (t).

In the example here addressed, 6 couples must be considered for reaching a linear convergence of 5 n Γ 32.

Figures 4 depict R 1 (x), R 2 (x), R 3 (x), R 4 (x) and figures 5 show S 1 (t), S 2 (t), S 3 (t), S 4 (t) for n Γ = 5 and n τ = 32. Figure 3 depicts the field u (x,t) the error u (x,t)u as (x,t).

Conclusion

The proposed approach transform an incremental BEM procedure into a decoupled one that needs the solution of some steady problems defined in space (Poisson equation in the case here addressed), and some ordinary differential equations that only involve the time coordinate. Significant reduction of CPU time is expected due to the non-incremental nature of the proposed technique, as well as

  as a function of S

Fig. 1 :

 1 Fig. 1: Colocation nodes (•), MLS nodes ( ) and integration points (+).

4 Table 1 :Fig. 2 :

 412 Fig. 2: Evolution of the error e n as a function of the space discretization and the number of functional couples, n RS , involved in the approximation.

Fig. 3 : 5 Fig. 4 : 5 Fig. 5 :

 35455 Fig. 3: u (x,t) et u (x,t)u as (x,t) for t = 0.2s and n Γ = 5


Patricio De La Cuadra

Benoit Fabre
email: fabreb@ccr.jussieu.fr

Nathalie Henrich

Thomas Robin

Analysis of jet instability in flute-like instruments by means of image processing : effect of the channel geometry on the jet instability

We are studying experimentally the behavior of jets in flute-like instruments. At the time being most models for jet oscillations are based on linear jet instability analysis described by Rayleigh. But there is no clear model for the initial perturbations induced by transverse acoustic velocity. Aiming to understand and describe the origin of the perturbation, we are currently investigating the jet motion under transverse acoustic perturbation generated by loudspeakers. The jet behavior is analyzed through flow visualization and image processing using a crosscorrelation technique. The paper will present results for different flue channel geometries. The linear instability is investigated. Experimental results are compared to a theoretical linear model.

Introduction

Instruments from the flute family share a principle of operation that consists of an unstable jet issuing from either a channel or the lips of a player. The jet flows in the direction of a sharp edge that we call the labium. The acoustic field due to the presence of the resonator triggers the jet instability and therefore the jet oscillates at the frequency of the acoustic field. The interaction of the perturbed flow with the labium provides the necessary acoustic energy to sustain the acoustic oscillation in the pipe.

Flutes present an immense variety of geometries and playing techniques. Among them, those that present a fixed geometry have particularly interested researchers in the last decades, mostly because of experimental simplicity.

The work presented here was triggered by the differences observed in the flue exit of various instruments. Organ pipes may show a flue channel ended by a 90 angle geometry. Recorders show a 45 cut at the flue exit and recorder makers insist on not rounding the edges of the flue exit, claiming that this would kill the instrument. This may seem contradictory since lip-played instruments obviously show a rounded flue exit.

The present paper explores the influence of the geometry of the flue channel exit (flue exit) on the oscillating behavior of a jet. Analysis are carried out using a nonintrusive technique based on image analysis of flow visualization obtained by the Schlieren method [START_REF] Merzkirch | Flow Visualization[END_REF]. We concentrate in this paper on the estimation of the convection velocity of perturbations on the jet for different channel geometries.

Jet formation and oscillation

The analysis of the oscillating jet is normally split into three parts: perturbation of the jet by the acoustic field (receptivity), propagation of the perturbation along the jet (instability) and generation of acoustic energy by the perturbed jet (sources).

Most receptivity models (e.g. [START_REF] Fabre | Physical modeling of flue instruments: a review of lumped models[END_REF], [START_REF] Fletcher | Wave propagation on a perturbed jet[END_REF]) are developed with the assumption of a potential flow where vorticity is conserved. Vorticity is generated in the boundary layers of the channel and injected into the flow at the flue exit. Perturbations are assumed to act locally at the point of flow separation, modulating the vorticity in the shear layers of the jet. None of the current models provide a description good enough to take into account the effect of the channel geometry (length and exit shape) on the jet perturbation.

The central speed of the jet U b at the channel exit can be estimated using Bernoulli's equation:

U b = 2p f ρ 0 (1)
Where p f is the pressure in the cavity before the channel, and ρ 0 the air density. Both the length of the channel and the flue exit geometry have an important influence in the velocity profile [START_REF] Segoufin | Experimental study of the influence of the mouth geometry on sound production in a recorder-like instrument: Windway length and Chamfers[END_REF]. Studying the frequency of edge-tone as well as visualizing transient attacks [START_REF] Segoufin | Experimental investigation of the flue channel geometry influence on edge-tone oscillations[END_REF], it has been observed that the convection velocity is higher on jets issuing from a short channel, corresponding to jets with sharper velocity profiles.

Due to the intrinsic instability of the jet, the oscillation is amplified while being convected downstream. And it travels at a velocity about one half of the jet speed [START_REF] Rayleigh | The theory of sound[END_REF]. In a linear model, we consider an amplification described by:

ψ = ϕ(y)e i(ωt-αx) (2)
Where ϕ is the complex amplitude of the stream function ψ, ω the angular pulsation, t time and α the wave number. Spatial analysis has been shown to provide a better fit to measurements, and can be written as a real ω, while looking for a complex wave number α = α r + iα i . The velocity of perturbation corresponds to c p = ω/α r while the spatial amplification is given by the coefficient α i .

Measurements of oscillating jets with short channel and squared exit have already been reported by Nolle [START_REF] Nolle | Sinuous instability of a planar air jet: Propagation parameters and acoustic excitation[END_REF] using a hot-wire. One of the benefits from choosing a visualization technique, as in the present paper, is that it is non-intrusive.

Experimental setup

Four flue exits are studied, of which the three with long channel are shown in Figure 1 (a fourth with squared exit and short channel is also included). Sequential images of the jet are taken with a digital camera [PCO,Sensicam, fast shutter] with exposure time set to one microsecond. Since the frequencies of jet oscillation are much higher that those provided by any modern camera, and because we want to be able to have an arbitrary number of images per cycle, a stroboscopic light is used. Its frequency is set such that multiples of its frequency are slightly out of phase with the excitation frequency providing an aliased representation of the oscillation. By changing this phase we can produce as many images per cycle as we want. Certainly they come from different cycles of the jet, but they can be collapsed to the same cycle since the oscillation is stationary.

Frequency of the camera was set to 14 fps while the acoustic excitation frequency was ranging from 70 to 1400 Hz to cover a Strouhal range from approximately 0.05 to 1.15 (Strouhal = f h/U b , with h the channel height and U b the jet speed at the exit deduced using eq.1). Approximately 100 images are taken covering two cycles of jet oscillation. Images are captured in raw bmp files, with size 1280 x 448, and 8 intensity bits.

The speed of the jet is controlled with a manual valve measuring the pressure in the cavity p f just before the jet formation. A Reynolds number based on the channel height h of 500 was chosen (Reynolds = U b h/ν, ν: kinematic viscosity) giving a jet velocity = 7.5 m/s which corresponds to normal blowing conditions for a recorder and assures a laminar behavior.

The amplitude of the acoustic excitation is measured through the acoustic velocity amplitude, which is created with two speakers in opposed phase, providing a more uniform field. A fixed amplitude was chosen for the whole experiment, and it was set to a value of 0.5% of the jet velocity. This is considerably lower that the 10% observed in real flutes but allow us to cover an interesting range of frequencies without over exciting the speakers.

The phase and amplitude of the acoustic excitation are measured using a velocity sensor (Microflown Technologies) placed as shown in Fig. 2. Due to its non-flat frequency response, phase and amplitude are compensated.

Image processing/data analysis 4.1. Cross-correlation method

Our goal is to determine the transverse displacement of the center of the jet as function of the distance from the flue exit for every image (Fig. 3). One image is defined as a reference, normally the first of the sequence. Then, the position of all other images are calculated as distances from that reference, obtaining a relative position for the whole sequence. Finally the position of the reference is determined and subtracted from the relative position to obtain the actual position. Cross-correlations are taken for the intensity values from the k th column of reference image X k ref , with that of the k th column of all other images from the sequence

X k t . Y k (l, t) = N i=1 X k ref (i) * X k t (l + i), l = -N + 1, . . . , 0, . . . , N -1 k = 1, . . . , P t = 1, . . . , M (3)
where, N is the number of pixels per column, P the number of columns per image, and M the number of images per sequence.

Then the position of the peak of Y k (l, t) is found, and the distance from the center calculated, giving the relative displacement. In order to estimate the jet deflection with a resolution better than 1 pixel, a parabolic interpolation is used, fitting y(x) = a * (x -p) 2 + b to the three highest points, and obtaining p as the interpolated value.

This method uses all the variability of the image sequence as valuable information to find the position of the jet. Fig. 4 illustrates the results obtained with it.

Processing data

All columns of the image oscillates in time at the frequency of the excitation. Therefore, for each column we fit a sinusoid at that frequency:

Y k (f ) = 1 N N -1 i=0 X k i e -j2πf f s (4)
The amplitudes |Y k (f )| and phases ∠Y k (f ) are shown in Fig 5 . The linear model described by eq.2, predicts exponential growth described by α i and a constant

Results and Discussion

The convection velocity of the perturbation c p for the four cases and complete range of frequencies is presented in Fig. 6, as fraction of the center jet velocity at the flue exit U b . In the case of the short channel with squared exit, the convection velocity is clearly faster than in the other cases, confirming the observations done by Ségoufin [START_REF] Segoufin | Experimental investigation of the flue channel geometry influence on edge-tone oscillations[END_REF] in the context of edge-tones. This can be interpreted, as suggested by Ségoufin, in terms of the boundary layers' thickness of the jet, which becomes thinner in the case of short channel. The magnitudes of the slowest three exits are in agreement both with the experimental measurements by Nolle [START_REF] Nolle | Sinuous instability of a planar air jet: Propagation parameters and acoustic excitation[END_REF] and theoretical values by Mattingly and Criminale [START_REF] Mattingly | Disturbance characteristics in a plane jet[END_REF].

For the three long channels, the convection velocity seems to be independant of the flue exit which confirms that it depends mostly on shear layer formation in the channel. A detailed analysis of the velocity field in the vicinity of the flue exit should be carried out in order to interpret these results further. The amplification coefficient α is presented dimensionless in Fig. 7. For each geometry exit, we observe an overall curved shape with a peak or preferred frequency that lies around Str = 0.4, which is in agreement with the theory [START_REF] Mattingly | Disturbance characteristics in a plane jet[END_REF] For the moment we have no tools to determine the position of the transition from exponential growth into a discrete vortex model [START_REF] Holger | Fluid mechanics of the edge-tone[END_REF]. Only the part of the jet that oscillates linearly has to be included in our calculations. The estimation of amplification α is very sensitive to the position of that boundary, producing data with bigger dispersion than the one from the propagation speed c p . Further work is needed in order to estimate the phase and amplitude of the excitation at the origin. An accurate model would need a detailed description of the total flow (from jet and from acoustics) at the flue exit.

Visualization of the oscillating jet together with the data analysis presented, provide a promising tool to explore some of the mysteries of the operation of flute-like instruments.

Figure 1 :

 1 Figure 1: Geometry of the three flue exits used, from left to right: Squared, Rounded, 45 degrees.

Figure 2 :

 2 Figure 2: Experimental setup. Jet exit, speakers and acoustic velocity sensor.

Figure 3 :

 3 Figure 3: Jet issued from a squared exit, Re: 500, Jet speed: 7.5 m/s, Excitation freq.: 532 Hz, Str: 0.44. Visualized with Schlieren technique: the two shear layers appear with opposite light refractive gradient, so that the upper shear layer gets darker while the lower shear layer gets lighter than the average gray level. The jet oscillating movement shows a very small transverse deviation at the beginning of the jet, slowly growing, up the point where the flow rolls up and finally breaks into vortices and later turbulence.

Figure 4 :

 4 Figure 4: Position of jet found with cross-correlation method, excitation frequency = 532.28 Hz, flue exit = squared. Time axis corresponds to frame times in seconds.

Figure 5 :

 5 Figure 5: Up: (thick) Amplitude of the detected sinusoids, (thin) Exponential fit. Down: (thick) Phase of detected sinusoids, (thin) Linear fit. Excitation frequency = 532.28 Hz, flue exit = squared

Figure 6 :

 6 Figure 6: Dimensionless velocity, flue exit short (dashed), squared(dashdot), rounded (solid) and 45 degrees (dotted).

Figure 7 :

 7 Figure 7: Dimensionless amplitude, flue exit short (dashed), squared(dashdot), rounded (solid) and 45 degrees (dotted).