Numerical investigation of vibration and dynamic pressure of a vertical axis wind turbine

In the environmental field, the problems of noise reduction have become a major preoccupation, particularly on the noise generated by the acoustic radiation pressure produced by wind turbines. This paper is aimed at presenting the investigation on the application of variational indirect boundary element method for study the acoustic radiation pressure produced by vertical-axis wind turbine. For this initiative, we considered Neumann boundary condition. The formulation has two advantages: the first one is to avoid the meshing of the fluid domain; the second advantage is to treat the singular integral of the Green`s function, solution of fundamental solution of the wave equation in frequency domain.


INTRODUCTION
In the environmental field, the noise pollution has become a major concern, particularly on the noise generated by the acoustic radiation pressure from wind turbines.In this sector, the audible noise is caused by the wind gliding over the blades and by the vibration of the generator (the noise reflects lost of energy).Modern wind turbines produce significantly less noise than older designs.Among these modern turbines, the vertical axis wind turbine is considered as a good solution to overcome the noise nuisance in urban environments where they are considered safer due to their lower rotational speeds.Vertical axis wind turbines can catch the wind from all directions at a lower wind speed than their horizontal axis counterparts.Indeed, vertical axis wind turbines are designed to spin at a slower speed than Horizontal Axis Turbines.Spinning at slower speeds allows the turbines to function at much higher wind speeds.It also reduces noise and vibrations, making Vertical Axis Wind Turbines a good alternative for noise reduction [1,2].
Most noise sources that are of concern to engineers can be modeled in terms of simple sources such as spheres, pistons in infinite bafflers, cylinders or combinations of these.In the case of wind turbines, the acoustic radiation pressure is modeled by experimental tests (acoustic and vibration).These techniques have remained empirical in that the design of the blades are determined largely by trial and error methods which are very costly as a lot of experience is required in order to estimate acoustic pressure radiation.It is therefore worthwhile to develop more efficient methods based, for example, on reliable numerical predictive tools to analyze and optimize the acoustic radiation emitted by wind turbines.This is enabled by the current development of precise sophisticated models that describe acoustic radiation and diffraction phenomena.
To accurately simulate the acoustic radiation pressure produced by a vertical-axis wind turbine, computational fluid dynamics in a moving fluid domain as well as fluid structure wind turbine interaction algorithms need to be performed.
In the present work, and as part of linear acoustics, we are interested in analyzing the acoustic radiation pressure produced by a vertical-axis wind turbine.The 3D variational indirect boundary element method and the Neumann boundary condition are considered for this purpose [3][4][5][6].The problem of interaction between a turbine structure (is assumed to be perfectly rigid) and the surrounding fluid is not considered.The fluid is supposed to be ideal and homogeneous medium.
The variational indirect boundary element method is extensively used in numerical acoustics, especially for solving radiation and scattering problems in infinite domains [3,7].This technique is based on integral equations that avoid the discretization of the fluid domain.Traditionally, methods based on integral formulation are known as boundary element methods.They generally use a collocation procedure [8][9] combined with a boundary finite element to generate full and asymmetrical system matrices, leading to inefficient treatment of the fluid and fluid-structure interaction problems [10][11][12][13].Moreover, these methods involve singular and hypersingular integrals (finite part) that are difficult to evaluate precisely thus leading to inaccurate results.The above difficulties can be circumvented by formulation of a variational boundary integral equation for the compressible fluid [3-4, 10, 14].
The main idea behind the formulation of a variational boundary integral equation consists of deriving a bilinear form defined on a bounded domain that leads to a positive definite quadratic form with a strong coerciveness property that guarantees the existence of a unique solution.However, this most desirable feature (strong coercivity) is not always attainable if the quadratic form satisfies the Garding's inequality.Fredholm's alternative may be applied so that uniqueness implies existence.The proposed formulation has several advantages: i) Over the finite element methods; it avoids the discretization of the fluid domain and ii) Over the collocation BEM formulation; it avoids explicit calculation of the finite part of singular integrals and leads after discretization by boundary element techniques to a small symmetrical algebraic system.
However, the BEMVI suffers from the singularity problem, which occurs when the double integral surface involves the same element.In previous years, effort have been made to calculate efficiently the singular integrals without sacrificing accuracy.In fact, special numerical treatment or analytical rearrangement of the resulting singular integrals is needed to eliminate difficulties related to their integrability.On the other hand, various approaches, mixing numerical and analytical quadrature methods have been successfully developed [3,4,11].In this work, we used the technique proposed by Alia et al, triangular element mesh [11].

GOVERNING EQUATIONS
For an harmonic disturbance of frequency "f " without any source or loss mechanisms, the pressure p satisfies the Helmholtz equation: where denotes the wave number, c is the sound velocity, w = 2pf is the pulsation and x is the field point position.For Neumann boundary condition which implies that the velocity is continuous across the boundary surface S: the pressure at any point within either one of the two acoustic domains can be expressed as following: ( where r is the fluid density, V n is the normal velocity, m = p 1p 2 is the double potential layer, n y is the unit normal at the location of the source point, is the Green's function with and is the distance between the field point x and the source point y.
For infinite domains we assume the pressure p satisfies the Sommerfeld condition: (4) Equation ( 3) states that the sound pressure at any point inside the acoustic domain can be obtained by integrating the equation on the boundary.This is the main idea of the BEM in which only the boundary information is needed to obtain the solution.However, equation ( 3) is not ready to be used because the double potential layer m on the boundary is unknown.In order to find it, we applied in what follows the boundary condition on the surface of the acoustic domain.By taking into account Neumann boundary condition given by equation ( 2), an integral equation can be derived for the velocity from equation (3).(5) This integral over the single surface S y , which is associated with the second derivative of the Green function, should be defined in the sense of Hadamard Finite Part (HFP) [4].
An increasing number of researchers have focused on hypersingularity and proposed both analytical and numerical techniques to handle it.A summary on regularisation methods for hypersingularity can be found in ref. [4].
In order to solve simultaneously interior and exterior problems, we will associate to equation (4) a variational formulation that presents a double advantage.First, it allows avoiding evaluation of HFP.In the other hand, it leads to a symmetric algebraic system.
The variational formulation can be derived by multiplying the integral equation [eq.5] by m(x) and integrating it over the BE model with velocity boundary condition on S.This formulation is based on the principle that the solution of the obtained equation will also minimize the following functional: The singularity which appears in the functional F can be reduced to a less singular form that is better suited to numerical calculations [1-3]: (7) By discretising equation ( 6) using BEM, the final numerical system is derived by imposing a stationary condition on F with respect to unknown primary variable m.Solving the obtained system allows to get the double layer potential at the boundary.Consequently, equation ( 3) can be used to compute the pressure at any point in the acoustic domain.

SOUND POWER
The sound power can be calculated by integrating the intensity over an imaginary surface surrounding the radiating object.Acoustic intensity, I n , is the time average of the rate of sound energy flow per unit area normal to the direction of propagation of the wave.It is a vector quantity in the direction of velocity.For time-harmonic waves, where the time dependence of pressure and velocity can berepresented by e -iwt , the intensity reduces to: (8) where * denotes the complex conjugate and Re indicates the real part.The radiated power W is given by following expression: (9)

DISCRETIZATION
In order to achieve a numerical solution of stationary functional (eqn.7), the surface is divided into finite boundary elements.The double potential layers μ(x) and μ(y) on the surface of the boundary element model are expressed as a product between their unknown nodal values and the element shape functions N. The discretized form of the functional J f may then be written as: (10) Imposing stationary condition on J f with respect to unknown primary variables μ: (11) The components of the elementary matrices A ij of the matrix [A] and the elementary vectors b j of the Vector{b} are given by the flowing expressions:  d d , ( , ) , square and symmetric since it is derived using a variational formulation.Once the double potential layer is calculated by solving the linear system given by Eqn (11), the acoustic pressure at any field point in the domain can be computed via integral equation (3).

NUMERICAL INTEGRATION
For the boundary finite element discretization of the fluid, we consider a triangulation T covering the boundary Γ, with linear triangular element L h , where h denotes the largest dimension of each element.For numerical integration purpose, four cases can be distinguished for a given pair of triangular elements , is the total number of boundary elements covering the surface Γ: Case 1: elements t i and t j are congruent Case 2: elements t i and t j have only one common edge Case 3: elements t i and t j have a common node (vertex) Case 4: elements t i and t j are disjoint.
Two types of integrals have to be evaluated during the formation of elementary vectors and matrices on the set of pairs of triangles For the case where elements t i and t j are disjoints and located at relatively large distance r from each other, we use standard Gauss-Radau integration formula [3][4]: (15) where m is the total number of integration points required to exactly integrate a polynomial of order on the triangle and is a regular function on its domain of integration corresponding to the product of kernels, shape functions and jacobian.
To evaluate the singular integrals, various approaches, mixing numerical and analytical quadrature methods have been successfully developed; useful references are cited in Reference [14].In this work, we proposed the Wang and Atalla [15] integration technique, proposed for quadratic elements and be extended to triangles elements by Alia and al [11].To demonstrate the accuracy and the efficiency of the numerical approach for evaluating the singular integral: (16) we consider there different situations: two triangles are coincident, two triangles have a common edge and two triangles have a common vertex.The results of the above and their situations are presented in the Table 1.From this table, when the integration is evaluated on the same element, the error is about 1.76%.However, it is about 4.26% when the singularity is concentered at the common edge of two adjacent elements and it becomes very weak (less than 0.2%) when the elements have one node in common.We can conclude that the presented method leads to good results.In what follows it will be seen that these errors do not affect the calculated acoustic pressure.

NUMERICAL VALIDATION
The indirect variational boundary element method (VBEMI) outlined in the previous section was implemented in the general purpose boundary finite element code developed by the author.This code was developed to study the acoustic radiation pressure produced by the three dimensional structure.All computations were performed on the PC in double precision.

ACOUSTIC RADIATION FROM A TRANSVERSELY OSCILLATING CIRCULAR DISK
For this application, we consider a circular disk of radius a = 1.0 m.The fluid medium surrounding the disk is air with sound speed c = 343.0m/s and means density r = 1.21 kg/m 3 .The disk is oscillating with a unit transverse velocity, perpendicular to the plane of the disk, in the normal direction; that is v n = -1m/s.The linear triangular mesh of the disk is shown in Figure 1.

426
Numerical investigation of vibration and dynamic pressure of a vertical axis wind turbine The analysis is performed for a frequency range of f n = {nc/2π} with n = {1, 2, 3, 4, 5}.The VBEMI model solves for the surface double potentials layer.These are used to compute surface pressures.
For validation, we consider the results in [16].In [16], the variational formulation derived for acoustic radiation and diffraction problems is solved using the Rayleigh-Ritz method.The basic functions are selected such that they satisfy certain characteristics of the exact solution for a disk.For example, the tangential derivative of the surface pressure is infinite at the edge of the disk and this knowledge is used in selecting appropriate basis functions.The treatment of this problem is described in details in Reference [16].
Figures 2.1 to 2.5 shows comparison of the surface real pressure and imaginary pressure from VBEMI method and from the results extracted from the numerical solutions published by Wu et al [16].The dimensionless surface pressures p/rv n c are plotted against dimensionless r = a for various ka = (1, 2, 3, 4, 5}, where r is the distance from the center of the disk.Figure 6.6 show the surface absolute pressure distribution from VBEMI.At the free edge, the pressures reduce to zero and the slopes of the curves tend to infinity.VBEMI solutions capture these characteristics well, and show good comparisons with the published numerical results [16].

PULSATING SPHERE
For the first case of validation, we consider the problem of a pulsating sphere of radius a = 1.0 m.The fluid medium surrounding the sphere is air with sound speed c = 343.0m/s and mean density r = 1.21 kg/m 3 .The wave number at a frequency ω is given as k = ω/c.The sphere is pulsating with a radial velocity v r = 1 m/s.The linear triangular mesh of the sphere is shown in Fig. 3.The pulsating sphere with a uniform radial velocity has exact solutions.For the exterior problem, the pressure at a distance r from the center of the sphere is given by: (29) Z is the characteristic impedance of air (Z = rc).For the interior problem, the pressure at a distance r from the center of the sphere is given by: (30)  The analytical expression for the radiated power by a pulsating sphere is given by the following expression: (31) The analysis is performed at a frequency f = 54.59Hzusing the VBEMI to solve for the double-layer potential m.This one, in turn, is used to compute field pressures at the specified points.The VBEMI solution on the sphere surface is (207.231± 0.029, -961.143± 0.115) for all nods, compares well with the analytical solution (207.515, -952.087).
Figure 4 shows field point pressure variation with distance from both VBEMI and Analytical methods.Since VBEMI solves both the interior and exterior problems simultaneously, we can compute field point pressures anywhere in the space from the double layer surface potentials.For r = a < 1 and r = a ≥ 1 the analytical solution from the interior and exterior problems is compared respectively with the VBEMI solutions.The comparisons show very good agreement between the solutions.The radiated power computed by VBEMI is 1292.91Watts, which matches well with the exact analytical solution of 1303.86 Watts for an exterior problem. Figure 5 depicts the variation of the radiated pressure with the radial distance computed analytically and by VIBEMI.Good agreement is observed between numerical and analytical solutions.
In a third example, let us calculate the radiated pressure for different frequencies at a point located at r = 4 m from the sphere center.From Figure 4, many peaks occur.They are related to the resonance of the interior volume of the sphere.These peaks (so-called critical or irregular frequencies) do not have any physical meaning at the considered frequencies.These special frequencies are related to eigen frequencies of the interior region and will be present if an interior BEMI problem is used.Since the interior region resonates and since both interior and exterior problems share the same integral operator, the integral equation governing the exterior problem should break down at the natural frequencies of the interior problem.Away from irregular frequencies, the presented numerical result shows good correlation with analytical solution as seen in Figure 6.The numerical algorithm for solving Numerical investigation of vibration and dynamic pressure of a vertical axis wind turbine irregular frequencies is not implemented yet in the code we developed to validate the presented the singular integration method.It is our goal to implement in the near future an algorithm to solve the irregular frequencies problem.The treatment of this problem is described in details in Reference [17].

ACOUSTIC RADIATION PRESSURE PRODUCED BY VERTICAL-AXIS WIND
This section is aimed at presenting the investigation on the application of VBEMI for study the acoustic radiation pressure produced by the experimental vertical-axis wind turbine of the Dermond Company * ; which is in an area belonging to the University of Quebec at Abitibi-Temiscamingue (see Figure 7).The geometric characteristics of the wind turbine are presented in Figure 7.The fluid medium surrounding the vertical-axis wind turbine structure is air.The vertical-axis wind turbine is excited at a frequency f by a normal velocity v n = V max (|n ⋅ i|).n is the unit normal vector, i is the unit vector in the x-direction and V max is the maximal amplitude of normal velocity.The remaining parts of the structure of wind turbine are assumed to be perfectly rigid.For numerical analysis, the wind turbine is replaced by a simple geometrical model and meshes using triangular elements (see Figure 7).In first step, for analysis, we will look the effect of the mesh element aspect on the results.For this, two different meshes are used: mesh 1 (2781 nodes and 4450 elements) and mesh 2 (2289 nodes and 3696 elements).The mesh1, involving 2781 nodes, is shown in Figure 8.In Figure 9, we presented the acoustic pressure radiated in exterior vertical-axis wind turbine obtained by the mesh 1 and mesh 2. According to this, which illustrates the results obtained, we note that both results coincides radiated pressure.This confirms the quality of mesh used.In the rest of the work, we consider only a mesh 1.
For analysis, four cases are considered to calculate the acoustic radiation pressure: 1.
On vertical plaque surface located at XZ plane (the surface plane is a square of side L = 5m); see Fig. 13.

EFFECT OF NORMAL VELOCITY AND FREQUENCY ON THE ACOUSTIC RADIATION
In this part we analyze two situations: the effect of the normal speed and frequency of the vertical axis wind turbine on the field of sound pressure (in dB).In the first case, the wind turbine is at frequency 1Hz.In the second case, the maximal normal velocity is (0.01 m/s).Figure 10 illustrate the effect of velocity (0.005 m/s and 0.01 m/s) of the radiated pressure and Figure 11 illustrates the effect of frequency (0.1 Hz, 0.25 Hz and 0.5 Hz).Note that the radiated pressure depending on the distance increases with velocity and frequency.However, the values remain low.

SOUND PRESSURE DISTRIBUTION BY VERTICAL-AXIS WIND TURBINE IN SPHERE SURFACE AND XZ PLANE
In Figures 12-13, we presented the acoustic pressure radiated in exterior of the wind turbine (sphere surface radiation distribution) and interior of the wind turbine wind in the XZ plane.We see that the noise generated on sphere-surface is very low for a distance of 20 m (the maximum is 3,3 Db).However, the acoustic pressure distribution in ZX plane is very important and the maximum is 87.9 db.

FAR-FIELD RADIATION PATTERNS OF THE VERTICAL-AXIS WIND TURBINE
Having solved the acoustic pressure on the surface of the wind turbine, we are able to predict the far-field radiation patterns.Here, we calculate the far-field acoustic pressure field from the rigid wind turbine at frequency 1 Hz.The values of the amplitudes of the far-field acoustic pressures given by VBMEI are plotted against the angle α which is measured from the central axis of the planes geometry projection wind turbine (XZ and YZ planes), such that α = 0 corresponding to the z axis direction (See figure 13).In order to depict the acoustic radiation beam patterns, we have plotted all values of the amplitudes of the far-field acoustic pressures in polar coordinates and normalized the result with respect to the maximal values evaluated at position (0, -6.366, 0.) in first far-field radiation patterns (YZ plane) and at position (6.761, 0., 0.) in second far-field radiation patterns in XZ plane.These results are presented at Figure 15.Numerical investigation of vibration and dynamic pressure of a vertical axis wind turbine Following this study, we found that the intensity of the noise structure is not significant.At this point, it will be useful in future work to induce rotation of the turbines in the analysis of acoustic radiation.Subsequently, a strong coupling fluid-structure vibroacoustic will be desirable.

(
12.a) (12.b)The global matrix [A] is resulting from assembling the elementary matrices A ij and the righthand side {b} results from assembling the elementary vectors b i .The complex matrix [A] is function on its domain of integration while is not necessarily regular.Utilizing local coordinates , integrals (24.a) and (24.b) can be written on the reference element as: (14.a) (14.b)

Figure 7 :
Figure 7: Geometric model of vertical-axis wind turbine

Figure 10 :Figure 11 :
Figure 10: Effect of normal velocity to sound pressure radiated by a Vertical-axis wind turbine (1Hz)

Figure 12 :Figure 13 :
Figure 12: Sound pressure distribution on sphere surface from a Vertical wind turbine (f = 1Hz)

Figure 15 :
Figure 15: Field point pressure comparisons

Figure 14 :
Figure 14: Plans YZ and XZ of wind turbine projection