A NUMERICAL INVESTIGATION ON THE VIBRATION OF A TWO-DECK EULER–BERNOULLI BEAM FLOODED BY A POTENTIAL FLOW

Summary This work concerns the structural vibration of a bladeless wind turbine, modelled by a two-deck Euler–Bernoulli beam, due to a surrounding potential flow. The deflection is governed by the Euler–Bernoulli equation which is studied first by a linear theory and then computed numerically by a finite difference method in space with a collocation method over the arc length, and an implicit Euler method in time. The fluid motion in the presence of gravity is governed by the full Euler equations and solved by the time-dependent conformal mapping technique together with a pseudo-spectral method. Numerical experiments of excitation by a moving disturbance on the fluid surface with/without a stochastic noise are carried out. The random process involved in generating the noise on the water surface is driven by a Wiener Process. A Monte Carlo method is used for stochastic computations. The generated surface waves impinge on the beam causing structural vibration which is presented and discussed in detail. By elementary statistical analysis, the structural response subject to the stochastic hydrodynamic disturbance caused by white noise is found to be Gaussian.


Introduction
A wind turbine is an engineering masterpiece that converts kinetic energy from the wind into electrical energy.It is eco-friendly with zero emissions and sustainable as renewable energy is collected.Many different types have been designed to improve the efficiency and reduce the manufacturing cost.The most common configuration consists of a large beam with blades atop.The reader may refer to (1-3) for a review.More recently, a novel project regarding vortexinduced bladeless wind turbines (BWT) was introduced by a tech firm called Vortex Bladeless.A comprehensive work by (4) studied the dynamic modelling of such context.In general, a BWT consists of a mast body mounted on a base via a flexible rod whose bottom is anchored to the ground.A schematic of a cylindrical BWT is presented in Fig. 1.In the article, we investigate the interaction of a BWT surrounded by a potential flow of shallow depth, which has potential applications in possible offshore installation.
The BWT can be considered as a classic Euler-Bernoulli beam whose history can be traced back to 1750 when two famous mathematicians (see (5)) introduced a linear partial differential equation (PDE), in one-dimensional space and time, which was used to characterise small deflections of a beam structure under the effect of a lateral load.Many different application areas, for example the structural design and analysis of cable-stayed bridges, roofs of football stadiums, high-rise buildings, and soil retaining walls, amongst others, with varying scales of loading.The motivation behind this work originated from the potential risk of structural damages caused by external strong impact forces such as those generated by floods, earthquakes or wind gusts.A good understanding of such vibration using mathematical analysis and computations has always been essential and greatly beneficial to the civil engineering community.Much work has already been achieved so far in the significant interest of wind loads (see for example (6)(7)(8) and the references therein).In this work, a valid numerical scheme based on a finite difference method by using the arc length as the canonical variable is introduced to compute extreme solutions with overhanging structures.Another important application is highlighted by beam-like structures in contact with water, for example, dams, longspan bridges, breakwaters and navigation locks, dam-reservoir systems and aircraft wings.It has been intensively studied in offshore engineering to understand the dynamical behaviour of the fluidstructure interaction, for example in the event of an earthquake, due to the interest of safe design and life time prediction.A comprehensive literature review can be found in (9).Early work described in the pioneering work of (10) in which the hydrodynamic loads were taken as an added mass attached to the interface between the structure and the water.Such a simplified model has been widely adopted in some codes of practice for example (11,12).Other work incorporated more complex features together with other numerical approaches such as a finite element or boundary element method (see for example (13)).A structural analysis was conducted in (14) with different settings of the surrounding water domain.In recent years, there has been an increase in research interest on this topic as there have been water-related catastrophes due to extreme weather leading to heavy rainstorms or tsunamis which in turn cause flooding and damages in low-level areas.To this end, the current work is focused on the flood-related problem in which the hydrodynamic loads are to be formulated by the full Euler equations.
Water wave has been a major research topic to scientists for centuries due to the common presence on the planet in the environment and daily life with multiple-scale features varying in the order of millimetres to the order of kilometres.To reduce the complexity of this flooded problem, a surface gravity wave is considered as the source of external load to the beam-like structure.The dynamics of the fluid surface are formulated under the framework of a time-dependent conformal mapping technique which was first pioneered in (15) and adopted by (16) in the context of water waves.The flow structure beneath the surface, including the hydrodynamic pressure (to be coupled with the beam structure), can be derived by elementary harmonic analysis and computed efficiently via the fast Fourier transform as shown in (17) and successfully applied to various problems (18)(19)(20).More recently, Flammarion and Pelinovsky conducted a statistical analysis of random water wave fields under the framework of the Benjamin-Ono equation in (21).For the present problem, it is also reasonable to include the effect of a stochastic disturbance in the dynamic boundary condition on the free surface to address the randomness in the surrounding environment.Stochastic computations are achieved in Monte Carlo simulations and followed by statistical analysis.
The purpose of the article is twofold: a. to model the beam-fluid interaction and conduct computations by a valid numerical scheme that deals with such coupling; 2. to understand the stochastic behaviour of the beam structure caused by random disturbances on the water surface.The article is structured as follows.The detailed formulation and numerical schemes are shown in sections 2 and 3.The computational results are presented and discussed in section 4. A concluding remark is made in section 5.

Formulation
A cylindrical BWT of height H is considered in a two-dimensional Cartesian coordinate system.It stands on flat ground and is surrounded by a potential flow of water with density ρ and mean depth h on both sides of the beam.Mean depth h is assumed to be small in comparison to H.The gravitational force acts in the negative y direction.At the fixed end on the bottom, the vibration of the base due to the pressure exerted by the water flow is no longer negligible as in the literature.Two different materials are considered for the base and the mast body with the former being much more rigid than the latter.Such physical feature is quantified by flexural rigidity or bending stiffness.All the structural properties are assumed to be uniform and time-independent.The notations are listed in Table 1.It is noted that the base is much more rigid than the mast.An illustrating representation of the BWT with various variables is presented in Fig. 2. The water on the left-hand side remains  stationary and is irrelevant to the current study.This is because the base of the BWT is significantly more rigid than the mast, and the maximum water depth, h, does not exceed the depth of the base.Consequently, in this context, the interaction between the water and the base can be approximated as that with a solid wall.The ratio of h over H is assumed to be small and denoted by in which δ 1.As the base is much more rigid than the mast, we assume that

.2)
A further assumption on the zero thickness of the beam is made.Non-dimensionalisation is achieved by selecting l, l/g, ρg (2.3) as the reference length, time and force per unit length.Under such scaling, the length of the base is fixed to be 1.All the other notations remain unchanged in the rescaled system.

Water waves
The associated velocity potential and streamfunction are denoted by φ(x, y, t) and ψ(x, y, t) respectively.The interface between the fluid and air is denoted by ζ (x, t).The fluid motion can be described by the velocity potential which follows ) ) All the other governing equations for the fluid remain unchanged.In the fluid body, the pressure field can be computed by the Bernoulli equation (2.8), which yields (2.12)

Euler-Bernoulli equation and the linear theory
The deflection w(y, t) of the BWT is governed by the rescaled Euler-Bernoulli equation (5, 28, 30) where Q(y, t) is an external load acting on the beam, that is the force per unit length, subject to the standard boundary conditions for a cantilever The coefficient a(y) and R(y) are defined by It can be easily seen that γ 0 γ 1 = O(δ 2 ).By separating the variables, the solution of the Euler-Bernoulli equation can be written in a series of modal compositions as follows where α j are the time coefficients and φ j are the spatial eigenfunctions or the so-called natural modes.By introducing the inner product .,. associated with the eigenfunctions defined by and making use of the orthogonality condition, it can be deduced that where φ j can be found to be in the form of Besides the boundary conditions (2.14), φ j , φ j , φ j , φ j from (2.20) and (2.21) are also matched at y = 1 giving 4 extra conditions.Therefore, an 8 × 8 linear system for the coefficients b ij is obtained.
The associated matrix has zero determinant yielding a solvability condition which is used to solve for the ω j , usually referred to as the natural frequencies of the structure.In the limit when δ → 0, the asymptotic behaviour of the φ j near the base is written by in which the first two terms are found to be zero by imposing the boundary conditions at y = 0 from (2.14).Here, we restrict Q to be the hydrodynamic pressure exerted by water.The vibration caused by aerodynamic pressure above the water surface can be computed and added into the solution subject to a linear superposition but is not considered in this work so that Q is only non-zero in the submerged region (y < h).In addition, we suppose that Q is of O(1), then there exists M ∈ R such that sup y<h |Q| = M, which can be used together with (2.1), (2.19) and (2.22) to derive As the deflection of the base is much smaller than the lengthscale of the water wave, it is reasonable to assume that the base part acts as a vertical wall whenever it interacts with external forces due to water flows.Also, the leading order term in the solvability condition for the natural frequencies is written by which is the characteristic equation of the natural frequencies for a uniform beam of unit length.
The ω j , up to the first-order accuracy, can be found by solving (2.24) using Newton's method.It is a straightforward exercise to derive the normal modes from (2.20) and (2.21), which are fundamental in structural analysis and engineering assessment.Next, we follow to investigate the dynamical response of the two-deck beam to the hydrodynamic pressure exerted by water waves for various values of water depth.It will be achieved by a numerical scheme coupling the beam deflection and the wave motion that is introduced in the following section.

Water dynamics
To compute the dynamics of the system (2.4)-(2.8),we employ a time-dependent conformal mapping, first pioneered by (15) and then adapted in water wave problem by (16), which conformally maps the physical domain into a simple geometry in the canonical plane.It has been successfully applied to many different problems concerned with the dynamics of a free surface in (23) for gravity waves, (22) for capillary-gravity waves, (24) for hydroelastic waves and (25) for electrohydrodynamic waves.In particular, the numerical experiments of excitation of gravity solitary waves were achieved numerically in (26,27).Meanwhile, the approach can be used to examine the flow structure in the fluid body given a free surface as shown in (18).
The pressure field in the fluid body can be then calculated from a given free surface y = ζ (x, t).The core idea behind such a technique is to transform the complex boundary in the physical domain [0, λ], where λ is the wavelength, onto a simple geometry in a new ξ − η plane that is a strip along ξ in [0, L] of finite depth d in our case.We denote the variables on the surface by The explicit form of the conformal mapping can be obtained by solving the boundary value problem in which = ∂ ξξ + ∂ ηη is the Laplacian in the canonical plane.Using elementary analysis, one obtains ) ) where .and T are defined by Here, 'PV' stands for Principal Value.The functions in the physical domain can be evaluated by ) and those in the canonical plane can be computed from the surface variables as follows ) (3.12) where F represents the fast Fourier transform.In terms of the new variables, the time-evolution equations become where J = X 2 ξ + Y 2 ξ is the Jacobian of the conformal map.The hydrodynamic pressure in the fluid body can be captured by using the Bernoulli equation (2.12) at an arbitrary level beneath the water surface.The reader may refer to (18) for more details regarding the harmonic analysis.
A fourth-order Runge-Kutta method is employed to compute the fluid dynamics via (3.14)-(3.15).For travelling waves with speed c, the system can be further simplified into one single equation The surface elevation can then be expressed as a Fourier series which is truncated after N terms ) are evaluated on these mesh points.The derivatives and the T-transform are computed via Fourier multipliers.The resulting system is solved by Newton's method.The l ∞ norm of the residual error is set to be less than 10 −10 for convergence test.Such a fast and accurate numerical method has been proven to be robust and efficient.

Beam deflection
The rescaled Euler-Bernoulli beam equation (2.13) is used to compute the beam deflection.The problem is parametrised over the arc length in the range [0, H].The N collocation points are chosen to be uniformly distributed with a step size H N along the domain The associated deflection is denoted by w i , that is w i = w(s i ).The boundary conditions in (2.14) can be imposed by introducing four ghost points as follows and setting The dynamics of the beam equation (2.13) is computed by an implicit Euler method in time and a second-order central difference scheme in space.The second-order equation (2.13) can be expressed as two coupled differential equations of the first order by using the artificial variable v = ∂w/∂t as follows where The computation starts at t = 0 and ends at t = T when the beam completes interacting with the incident water wave.In practice, T is selected to be 200, and the temporal step size t is chosen to be T n.The discretised variables are denoted by

.25)
The matrix form of the backward Euler scheme can be written in the form of where I N is the identity matrix and M is the matrix containing finite difference replacement of the fourth spatial derivative, characterised by central difference coefficients: 6 on the diagonal, −4 on the adjacent diagonals, and 1 on the second diagonals, all scaled by the factor 1/(H/N) 4 .The two vectors T are defined similarly as in w i .Multiplying both sides by the inverse of the matrix from the left-hand side of (3.26) completes one temporal increment.The numerical stability is guaranteed due to the nature of the implicit method.
The efficient performance of this scheme and its advantages over the finite element method were discussed in detail in (28).

Summary
To summarise, a fast hybrid scheme is introduced to solve the problem of fluid-structure interaction of a beam with its base immersed in an inviscid incompressible flow.The fluid dynamics is computed by a Fourier spectral method combined with a fourth-order Runge-Kutta method whereas the beam deflection is simulated by a finite difference scheme with collocation over the arc length together with an Implicit Euler method.The two parts are coupled by the load term Q that is obtained by by considering the pressure difference exerted on the two sides of the beam structure.The pressure exerted on the left side of the beam is constant since the water level on that side remains unchanged under the assumption that the beam acts like a vertical wall during the beam-fluid interaction.The pressure beneath the surface exerted on the right side of the beam can be computed from the Bernoulli equation given a surface profile as illustrated in (18).The elastic collision between the fluid and the non-deformable submerged beam is achieved by using the method of images in an extended domain [0, 2L] where two identical wave dynamics are simulated in [0, L] and in [L, 2L] respectively, propagating in opposite direction towards the two boundaries x = 0 and x = 2L.The elastic collision takes place at x = 0, where the two wave dynamics meet due to the periodic boundary from the Fourier spectral method.In all the computations, the step sizes are selected as ξ = 0.1, s = 0.02 and t = 0.01 for achieving high accuracy.The parameters of conventional materials for a wind turbine are as follows which will be fixed in the subsequent computations.

Dynamics
Numerical experiments conducted on excitation are described in this section.The free surface is initially flat and under the effect of a single moving disturbances which is placed on the right hand side of the beam at a distance and is switched on at the launch and switched off at a given time later.The excited waves continue to propagate towards the beam and are reflected due to a collision.The computations are stopped long after the end of fluid-beam interaction.To form nonlinear solitary waves, fully localised forcing such as a normal distribution (see ( For the integral (4.1) and (4.2), the integration interval from 0 to H represents the entire length of the beam.The integral sums up the absolute values of the curvature over the entire length of the beam, providing a measure of the total bending at time t.The absolute value is used to ensure that all displacement and curvature contributions are positive, emphasizing the magnitude of the bending without regard to the direction of the displacement and curvature.
4.1.1Single pressure distribution.A single pressure distribution is considered as the disturbance on the water surface in this subsection.It replaces p = 0 in (2.8) on the water surface and is initially placed at x = d moving leftwards with a speed U under the form of It is switched on at t = 0 and later off at t = 20.Taking out the pressure at a certain point generates a clean solitary wave rather than a periodic shedding of solitary waves (see (26)).The value for d is chosen to be 80, far away from the beam-like structure, such that there is enough space for the excited waves to propagate prior to the fluid-structure interaction.The wave dynamics after taking out the pressure disturbance are to be investigated carefully.The velocity U is selected to be 1.In the case where h = 0.5, the snapshots during the excitation stage are shown in the top two graphs from Fig. 3.As can be seen from graph (a), a stable depression solitary wave is excited.The excited waves continue to move towards and then collide with the beam situated at x = 0 as in graph (c) of Fig. 3.After the interaction, the waves bounce back and travel in the opposite direction as depicted in graph (d).Qualitatively similar surface displacements due to the moving disturbance have been observed from the numerical experiments with the same setting for other water depths such as h = 0.2 and h = 1.We omit the details here.The associated beam deflections for several time domains, or simply the range of the beam deflections, are depicted in Fig. 4. The graphs exhibit very similar features.To make further comparisons of the structural behaviours, the time-evolution of w † and κ † are sketched in Fig. 5.These two quantities remain zero until the propagation of the excited water waves reaches the beam structure.The first interaction times obtained from the three experiments for h = 0.2, 0.5 and 1 are found to be identical.Also, w † and κ † both exhibit oscillatory behaviours at a frequency that has a weak dependency on the value of h.This is a piece of strong numerical evidence from Figs 4 and 5 that the responses from the beam structure due to the interaction with the water wave are qualitatively similar for various water depths, and only the magnitude is subject to a factor depending on the value of h.In the subsequent computations, the water depth is fixed to be 0.5.

Stochastic disturbance.
A stochastic noise is included in the moving disturbance (4.4) as follows  where χ (t) is a standard white noise with respect to time written in the form of a generalised time derivative of a Wiener process.By (29), it can be discretised as where N(0, 1) denotes a normally distributed random variable with zero mean and unit variance.It is noted that the property has been used in (4.6).The parameter σ from (4.5) controls the amplitude size of the extra random term or its standard deviation in a statistical sense.The stochastic effect is introduced in order to address the complex effect of a raindrop or wind force and therefore assumed to be much smaller than the deterministic pressure.The parameter σ is fixed to be 0.001 in all the stochastic simulations.The stochastic term is assumed to be uniform in space for the reason that the length scale of this problem is small in comparison to the meteorological scale.Again, the pressure disturbance moving at a unity speed is turned on at t = 0 and off at t = 20.A thousand simulations were performed in one stochastic example.To reduce the computational costs, the temporal domain is shortened to [0, 120].
The beam deflection purely due to the stochastic disturbance is measured by  by subtracting the deterministic response (w 0 and κ 0 ) from the computational outcome.The profiles of ŵ at different times are depicted in Fig. 6.The behaviours of δ 1,2 from 1000 simulations are demonstrated in Fig. 7.It is interesting to observe that the mean value of δ 1 and δ 2 oscillate and decay, resulting in the beam returning to its undisturbed state in the long term.Meanwhile, δ 1 and δ 2 are less than 5% of the deterministic response w † 0 and κ † 0 .The medians oscillate around and eventually approach zero as time runs on whereas the interquartile range generally decreases in time t.The outliers with an appearance on only several occasions are not a concern as they are still of small size at the order of 10 −4 .Overall, the computational results illustrate that the stochastic noise cannot cumulate enough energy to make a significant impact on the beam causing a highly overhanging structure or deflecting over the design limit that is usually 0.4% of the beam length.Increasing the value of σ could completely change the story as the stochastic disturbance becomes the dominant external forcing imposed on the water surface, which violates the assumption to capture the physical feature of wind load or raindrop in the problem of a flooded beam-like structure.

Statistical analysis
To provide a rigorous insight into the vibrational behaviour of the beam, statistical analysis is conducted.The primary objective is to ascertain the nature of these vibrations subject to stochastic hydrodynamic disturbance and understand their behaviour over time.Again, the statistical analysis is based on examining 1000 simulations subject to the stochastic hydrodynamic disturbance in which h = 0.5 as previously presented.
The initial phase of the analysis involved the computation of the mean and standard deviation of the beam vibration ŵ at the designated times for all the points.Figure 8 shows the deduced mean amplitude of the vibration of the beam is of the order 10 −8 .It indicates the inherent inclination of the beam to resist deviations from its equilibrium state with stochastic disturbance.The computed standard deviation is of the order 10 −7 .It suggests a relatively tight cluster of values around the mean, suggesting a consistent and stable response of the beam to disturbances.
The distribution of the beam vibration ŵ is revealed through the probability density function (PDF) in Fig. 9 at several points (with heights 2, 6 and 10 of the beam) of the beam at the end of the simulations.The PDF plots exhibit a normal distribution tendency.It implies that the stochastic beam response is governed by a plethora of small-magnitude random disturbances.Further validation of the Gaussian nature of the vibrational data is achieved through the quantile-quantile (Q-Q) plots.The linear trajectory observed in these plots emphasised the normal distribution characteristics   inherent in the beam vibration.The Gaussian nature of these vibrations could be attributed to the cumulative effect of numerous small independent forces acting on the beam, which is a characteristic inherent to the Central Limit Theorem.
Beyond the basic descriptive statistics, the higher-order analysis of the displacement distribution of the beam is assessed by investigating the skewness and kurtosis to provide a more holistic perspective on the vibration distribution as shown in the bottom panel of Fig. 8.The skewness values, ranging between (−0.25, 0.5), suggest an approximated symmetric distribution.A skewness value An important aspect of understanding the beam response is to examine the spatial correlation, specifically the covariance between the displacements of different points on the beam at a given time.The covariance shown in Fig. 10, of the order 10 −15 , suggests an extremely weak correlation.This weak spatial correlation signifies that the displacement at one point on the beam is largely independent of the displacement at another point.This is a key observation as it indicates a uniform and independent response to the external disturbances across the length of the beam.

Kolmogorov-Smirnov test.
To validate claims about the normal distribution of beam vibrations, we employ null-hypothesis testing with the one-sample Kolmogorov-Smirnov (KS) test for a 5% significance level, which examines whether a sample conforms to a specified theoretical distribution, for example the normal distribution.This test is practically recommended for sample sizes of n 50 to assess data set conformity with a normal distribution.The test statistic, denoted by d, measures the distance between the cumulative distribution function (CDF) of the data and the CDF of the reference distribution.The associated p-value quantifies the probability of obtaining the test result as extreme as the data actually observed under the assumption of the null hypothesis.A p-value exceeding 0.05 (the chosen significance level) is a common threshold to reasonably infer that the data follows a normal distribution.We perform the KS test directly on the data of the original displacement w.As the KS test also assumes that there are no tied values (identical or very similar data points) in the data set since tied values can compromise the accuracy of the test causing inaccurate outcomes, duplicate values are removed from the data set to mitigate this issue.This step of data cleaning is crucial for obtaining reliable test results.Another advantage of such a sampling technique is reducing the sample size from 1000 to approximately 100.
We follow to evaluate the p-values for the null-hypothesis significance test at several positions on the beam at the final time as shown in Fig. 11 for y = 4, y = 6 and y = 10.They are found to be greater than the threshold value, and so are at other heights on the beam.It exhibits the failure to reject the null hypothesis, which strongly suggests that the beam displacement w closely conforms to a normal distribution.

Conclusion
In this work, the problem of a beam flooded by a potential flow was considered.The dynamic behaviour of the beam colliding with water waves was computed numerically.The numerical scheme was successfully implemented to study the beam-fluid interaction.Several excitation experiments in the absence/presence of randomness were examined.The computational results were analysed and discussed in detail.The structural response of the beam to the stochastic disturbance governed  by the white noise was shown to be Gaussian by regression and hypothesis testing.Engineers and researchers can now predict the beam behaviour under certain conditions, allowing for better design strategies and mitigation techniques.Moreover, this statistical insight complements previous findings, offering a more holistic understanding of the beam vibration under the influence of stochastic disturbances.

Fig. 2
Fig. 2 Schematic of the problem.
.20) φ j = b 21 cos δω j y + b 22 sin δω j y + b 23 cosh δω j y + b 24 sinh δω j y , for 0 y < 1.(2.21) ) where the coefficients a n , b n are unknowns.In all the computations, N = 2048 Fourier modes are used.The domain [0, L] is uniformly discretised into N mesh points, and L = 200 is commonly used in the computations.The dynamic equations (3.14)-(3.15

Fig. 3
Fig. 3 Snapshots of the water surface at t = 15, 32, 80, 120 plotted in the physical plane in the excitation experiment under a pressure distribution defined in (4.4) with A = 0.03, d = 80 and U = 1 for water depth h = 0.5.The disturbance is turned off at t = 20.

Fig. 5
Fig. 5 Time-evolution of w † and κ † in the experiment under the effect of a single moving disturbance (4.4).

Fig. 6
Fig.6 Profiles of ŵ at specific times for water depth h = 0.5, plotted in the physical plane.Grey lines are the beam deflections for 1000 simulations at specific times; black lines are the profiles averaged over 1000 simulations.

Fig. 7
Fig. 7 The values of δ 1 (left) and δ 2 (right) versus t obtained from 1000 different simulations.The mean values of δ 1 and δ 2 are plotted in black lines.

Fig. 8
Fig.8 The profiles of the first four statistical moments at t = T for all the points of the beam.

Fig. 9 (
Fig. 9 (Top) Histograms and the fitted probability density function.(Bottom) The quantile-quantile plots drawn in grey dots (bottom panels) at several positions of the beam when t = T.The black lines are the reference lines with a unity slope.
Z. LIU ET AL.