On a Novel Approximate Solution to the Inhomogeneous Euler–Bernoulli Equation with an Application to Aeroelastics

: This paper focuses on the development of an explicit ﬁnite difference numerical method for approximating the solution of the inhomogeneous fourth-order Euler–Bernoulli beam bending equation with velocity-dependent damping and second moment of area, mass and elastic modulus distribution varying with distance along the beam. We verify the method by comparing its predictions with an exact analytical solution of the homogeneous equation, we use the generalised Richardson extrapolation to show that the method is grid convergent and we extend the application of the Lax– Richtmyer stability criteria to higher-order schemes to ensure that it is numerically stable. Finally, we present three sets of computational experiments. The ﬁrst set simulates the behaviour of the un-loaded beam and is validated against the analytic solution. The second set simulates the time-dependent dynamic behaviour of a damped beam of varying stiffness and mass distributions under arbitrary externally applied loading in an aeroelastic analysis setting by approximating the inhomogeneous equation using the ﬁnite difference method derived here. We compare the third set of simulations of the steady-state deﬂection with the results of static beam bending experiments conducted at Cranﬁeld University. Overall, we developed an accurate, stable and convergent numerical framework for solving the inhomogeneous Euler–Bernoulli equation over a wide range of boundary conditions. Aircraft manufacturers are starting to consider conﬁgurations with increased wing aspect ratios and reduced structural weight which lead to more slender and ﬂexible designs. Aeroelastic analysis now plays a central role in the design process. Efﬁcient computational tools for the prediction of the deformation of wings under external loads are in demand and this has motivated the work carried out in this paper.


Introduction
Pressures to minimise the environmental impact of air travel have led to aircraft manufacturers becoming focused on reducing the fuel consumption of the next generation of civil transport aircraft. Manufacturers are starting to consider configurations with increased wing aspect ratios and reduced structural weight, and this is leading to more structurally flexible designs [1][2][3][4]. Engineers can no longer assume that the airframe is rigid, and it becomes increasingly important to take account of the aeroelastic behaviour earlier in the design process. Efficient computational tools for use at the conceptual design stage for the prediction of the deformation of wing beams with representative physical properties under varying external loads are now in demand. Such tools are also applicable to the aeroelastic analysis of helicopter rotor blades and wind turbines. Recent research works related to aeroelastic analysis were conducted by Amoozgar et al. [5], Vazhayil Thomas et al. [6], Rajpal et al. [7], Tsushima et al. [8], Rong et al. [9], Sommerwerk et al. [10] and Qin et al. [11]. It is with these developments in mind that the present work was undertaken.
The Lax-Richtmyer stability criteria were applied to the scheme to derive an upper limit on ∆t to ensure numerical stability. The concern is that such a restriction on the size of ∆t might render this finite difference scheme computationally inefficient. However, this concern has proved unfounded in practice. To strengthen confidence in the validity of the scheme and its implementation, we present the results of a grid sensitivity study using the Richardson generalised extrapolation method (as can be seen in the work of Roache in [66]) and these results demonstrated that the scheme is indeed grid convergent.
We followed this by presenting the results of three sets of computer experiments. The first set simulated the behaviour of the unloaded beam. Careful choice of the initial and boundary conditions allowed us to compare the results with the homogeneous analytical solution provided in the earlier part of this paper. The second set of computer experiments simulated the time-domain behaviour of a beam with properties chosen to represent the slender flexible wing structure of a generic aircraft. This consisted of the stiffness and mass distributions varying along the length of the beam, while subjected to inertial and fluid-dynamic loading during a gust encounter. It is crucial to appreciate that in this work the external load distribution due to the fluid-dynamic lift is dependent on the surface geometry, which in turn depends on the beam deformation. In this way, the deformation feeds back into the external load distribution. The third set of computer experiments was used to validate the finite difference method against experimental data. Experiments were conducted as part of background research at Cranfield University to measure the static deformation of the tip of a horizontally oriented aluminium alloy beam cantilevered at one end as different transverse point loads are applied to the opposite end. Computer simulations were conducted with the replicated boundary conditions and the results were compared with the experimental results. This paper concludes by summarising the main findings and giving suggestions for further work to expand the applicability of the methods presented.

Governing Equations and Solution Methodology
In this section, we develop a numerical approach for solving the inhomogeneous Euler-Bernoulli equation. The derivation of an analytical solution to the homogeneous version of the Euler-Bernoulli equation with simple boundary and initial conditions is provided in Appendix A which we use in Section 3 to validate the numerical scheme developed in this work. For the inhomogeneous equation, when the function describing the transverse loading is taken into account, the unsteady numerical method was employed and a stability analysis was also carried out. We assumed that (a) deformations remain elastic which means that the deformations are small and that strain in the material remains within the elastic range; (b) plane sections remain plane, which means that cross-sections do not distort out of plane during bending; (c) the material properties are isotropic; and (d) the beam is not subject to any longitudinal loading. These assumptions restrict the validity of the Euler-Bernoulli theory to those cases where the deformation is small relative to the length of the beam. This means that nonlinearities associated with large deflections are neglected in the work presented here. In such cases, the beam bending theory of Timoshenko can be considered.
Under the assumptions stated above, the Euler-Bernoulli equation can be used to predict the dynamic deformation behaviour of a slender beam under a laterally applied distributed load. The inhomogeneous Euler-Bernoulli equation in a rectangular Cartesian coordinate system (see Figure 1) can be written as where E = E(x) is the modulus of elasticity of the beam material, I = I(x) is the second moment of area of cross-section of the beam taken about the neutral axis parallel to the plane z = 0, µ = µ(x) is the mass-per-unit length along the beam, the function w = w(x, t) is the lateral deflection of the beam as a function of distance along the beam x and time t, and q = q(x, t) is the given external force per metre applied to the beam as a function of distance along the beam x and time t.
The coordinate system of a structural beam along with the notation used in this paper.
In real-world engineering applications, mechanical damping causes the amplitude of oscillation of a beam under a steady externally applied load to decrease over time and the deformation tends towards an equilibrium deformation profile. To capture this behaviour, we included a velocity-dependent damping term in our simulation. We assumed that the mechanical damping force was linearly proportional to velocity, and directly the opposite to the local deflection velocity. Including the damping term and re-arranging, we have: where β is the non-dimensional damping coefficient.

Development of a Numerical Approach for Solving the Inhomogeneous Euler-Bernoulli Equation
In this section, we develop a numerical scheme for solving the inhomogeneous Euler-Bernoulli Equation (1) by using finite difference approximations. We emphasise that the proposed approach overcomes the assumption that the deformation function is the linear superposition of basis mode shapes, and thus we avoid the necessity of defining the shape of the mode basis functions at the outset. The numerical scheme is applied to a single beam experiencing a transverse load distribution q(x, t) in a single direction (see Figure 2). We note that the load distribution q(x, t) is not a follower distribution, and therefore the method described here is limited to conditions where values of ∂w/∂x are small.
Geometric nonlinearities were taken into account in this work by permitting arbitrary distributions of structural properties as a function of distance along the beam x, and these distributions are sampled at discrete positions corresponding to the beam nodes of the finite difference method described below. These variable structural properties are E = E(x); I = I(x); µ = µ(x); and q = q(t, x). If the distribution of these beam properties rapidly changes, our approach would be to increase the fineness of the discretisation and increase the number of nodes in the finite difference scheme. The purpose of the time-marching finite difference scheme is to determine w(x i , t n+1 ), given w(x i , t n ), and w(x i , t n−1 ) for all i ∈ {0, . . . , N}, where we have adopted the subscript n in keeping with the convention used in finite difference schemes for indexing the time levels. We combined finite difference approximations for the various derivatives in Equation (1) to form the finite difference version of the Euler-Bernoulli equation and we re-arranged the equation to obtain an expression for w(x i , t n+1 ) in terms of w(x i , t n ) and w(x i , t n−1 ) for all i ∈ {0, . . . , N}.
We emphasise that the function w(x, t) represents the deformation of the initially straight beam, and is measured perpendicularly to the straight beam. Since no account was taken of the change in the angle of the beam as w varies, this approach necessarily ignores the preservation of beam length. This is an additional reason for the imposition of the small ∂w/∂x condition of the method.
We assume that the beam is initially at rest so that: Using numerical techniques to propagate the geometry of the beam from one discrete time level to the next involves solving the Euler-Bernoulli equation for the unknown node displacement at each node. It is important to note that great care must be taken in computing the numerical approximation of the derivatives in Equation (1). To perform computations for nodes near the ends of the beam, it is necessary to fill the stencil beyond the ends with artificially generated values of w. We call these artificial values "ghost" points and we call the artificially generated nodes at which these values are taken "ghost" nodes. We will show that contrarily to being an impediment to the method, the values of w at the ghost points are fully determined by the boundary conditions.
We proceed to substitute for partial derivatives in the inhomogeneous Equation (2) using finite difference expressions. Without loss of generality, we assumed ∆t is a constant for all time levels and ∆x is a constant for each point along the beam. Evaluating the inhomogeneous Euler-Bernoulli Equation (2) at the discrete coordinates x i and t n , we obtain: Substituting central second-order accurate temporal and spatial discretisation schemes, we can write the finite difference approximation of Equation (4) as where x n+1 = x n + ∆x, etc., and where S is defined as where the superscripts b, c, and f denote backward, central, and forward second derivatives to third-order accuracy, respectively, and where F(x i ) = E(x i )I(x i ). For the three second-degree spatial derivatives in Equation (6), we employed third-order accurate finite difference discretisation schemes, and we write: where the central third-order accurate spatial second-degree derivative in Equation (6) is estimated by the backward third-order accurate scheme is estimated by and the forward third-order accurate discretisation is estimated by The process of forming the finite difference representation of a partial differential equation requires that each derivative in the original equation be replaced by a finite difference approximation of a given order of accuracy. The order of accuracy directly affects the number of terms in the Taylor expansion that are taken into account. It can be a laborious task to derive these approximations by hand. We used an algorithm suitable for computer automation in the derivations of the expressions in Equations (8)- (10). Such algorithms are frequently taught in university courses but we did not find them in the published literature. For this reason, we included a derivation of the algorithm in Appendix B. A Matlab ® script that implements this algorithm can be found in the MDPI code database (see Appendix E).
A specific combination of left, centre, and right finite differences are contained in Equation (6). This choice of biases is largely based on our computational experience. We tried to use other combinations of biases, but these attempts exhibited unexpected oscillatory characteristics that do not correspond to physically realistic behaviour found in the experiment. For example, replacing the biases in Equation (6) with centred finite differences throughout results in the wild unbounded continual oscillatory behaviour of the function w that appears right from the start of the simulation. The authors believe that the suitability of the choice of biases in the different terms is associated with the route of information propagation through the algorithm. To illustrate our heuristic justification, we present plots of stencils of the nodes used in two different simulations (see Figures 3 and 4). A derivative sampling scheme that gives rise to wild oscillatory behaviour is presented in Figure 3 and the derivative sampling finally chosen for this work is shown in Figure 4.  It is important to note that the oscillatory scheme (see Figure 3) required the use of third-order accuracy for the approximation of both second-degree derivatives, while the physically realistic scheme (see Figure 4) only requires third-order accuracy for the approximation of the first second-degree derivative. This suggests that taking higher-order accuracy alone does not lead in itself to the better capture of physical phenomena.
Re-arranging Equation (5) and dropping the order terms gives: where S(x i , t n ) is given by Equation (7). The stencil of Equation (11) is illustrated in Figure 5. The outputs of the algorithm are approximations of the deflection function w at discretely sampled intervals along the length of the beam x and at discretely sampled times t greater than 0. In the next subsection, we derive sufficient criteria for the finite difference scheme to be stable.

Extension of the Lax-Richtmyer Stability Criteria to the Fourth-Order Euler-Bernoulli Equation
There are several established methods for determining conditions on the spatial step size ∆x and the time-step size ∆t in PDEs such that instabilities are avoided. An existing method suitable for use in explicit parabolic finite difference schemes of second-order linear PDEs is given in [67].
In order to assure the stability of the fourth-order Euler-Bernoulli equation presented in this paper, we propose to extend the application of the Lax-Richtmyer stability criteria to higher-order schemes to obtain an upper limit on ∆t which is dependent on ∆x. We focus on the coefficient of the w(x i , t n ) term in Equation (11). According to the Lax-Richtmyer criteria, the finite difference scheme will likely be stable if the coefficient has a positive pre-sign.
First, considering the multiplicative term in Equation (11), we see that: for all β > 0 and all ∆t > 0, so it remains to consider the coefficients inside the square braces in Equation (11). Grouping all terms in w(x i , t n ) and explicitly writing this condition, we require that: which leads to the inequality for a stable time-step size of the discretised Euler-Bernoulli equation given by Since this condition must hold for the whole beam, we see that ∆t must satisfy: The Lax-Richtmyer stability criteria states that the pre-signs of all coefficients of w depending on ∆x and ∆t on the right-hand side of Equation (11) should be positive. Taking the coefficient of w(x i−4 , t n ) as an example, we have: which is always positive for all ∆x and ∆t by virtue of the fact that F(x j ) is positive for all j ∈ {1, . . . , N}. Similar expressions exist for the coefficients of w( , and w(x i+3 , t n ), the coefficients are negative irrespective of the choice of ∆x and ∆t, and so the Lax-Richtmyer stability criteria can never be satisfied by these terms. It may appear that this jeopardises the stability of the algorithm. However, some terms with a negative pre-sign are in fact necessary; otherwise, the function w(x i , t n ) would be a monotonic increasing function of x i and t n .

Specification of Boundary Conditions and Values of Functions F and W at the Ghost Nodes
As illustrated in Figure 6, to compute w(x i , t n+1 ) for i ∈ {1, . . . , N} in Equation (11), two ghost values of F are required beyond each end of the beam. We defined these values by linear extrapolation as and: For the cantilevered beam considered in this paper, the physical boundary conditions on w are: To ensure that the gradient of the w function vanishes at x 0 , we set w(x −1 , t n ) = w(x 0 , t n ) = w(x 1 , t n ) = 0 for all n ∈ N. We also set w(x −2 , t n ) = 0 for all n ∈ N so the first and second left hand ghost points are determined by the second boundary condition. Considering the third boundary condition, we have: and re-arranging this gives: so the first right ghost point is determined by the third boundary condition. Finally, we consider the fourth boundary condition, and we require: and re-arranging this gives: so the second right ghost point is determined by the fourth boundary condition. The third and fourth left and right ghost points remain to be determined. We arbitrarily constrain the third and fourth left-hand ghost point values of w to be zero. We set the third and fourth right-hand ghost points using linear extrapolation as and: The finite difference scheme described in this paper together with the boundary conditions were implemented in a Matlab ® script which is included in the dataset linked to this paper (see Appendix E). All simulations were performed using this script.

Grid Sensitivity Study for the Fourth-Order Euler-Bernoulli Equation
We begin by showing that the finite difference scheme for approximating the fourthorder Euler-Bernoulli equation with damping is grid convergent. We present the results of a grid sensitivity study using the Richardson generalised extrapolation method as detailed in [66]. For this activity, we used a uniformly applied constant external load function q(x, t) = q, a constant beam structural properties E(x) = E and I(x) = I, and a damping coefficient β = 10 throughout. We examined the accuracy of the finite element scheme as ∆x is varied. For each choice of ∆x, the value of ∆t is chosen such that the Lax-Richtmyer criteria were met for the w(x i , t n ) term for all i ∈ {0, . . . , N}. We performed four simulations of the finite difference scheme and we identified variables associated with each of the different simulations by the subscript j.
We selected four values of ∆x denoted by ∆x j for j ∈ {1, 2, 3, 4} such that: so that, in general, ∆x j+1 ≈ ∆x j /2 for each j ∈ {1, 2, 3}. The number of nodes for each j is given by and in order to keep the mass of the beam M constant and independent of j, the node mass M j (x i ) = M j is taken to be: Now, we select a fixed time t > t 0 where t 0 is the time at the start of a simulation, and a fixed C ∈ (0, 1). Then, we define a node index i j such that i j /N j ≈ C for each j ∈ {1, 2, 3, 4}. In this way, the i th j node is at approximately the same point along the beam for all values of j. We proceeded to perform the four simulations of the finite difference scheme, one for each ∆x j , and consider the predicted value of the function w at the point (x i j , t ) for each simulation j. For a more convenient notation, we set: for each simulation.
Since the Lax-Richtmyer stability inequality in (15) contains a (∆x) 2 multiplicative term, it is clear that, as the value of ∆x reduces, a smaller step size in t must be chosen to ensure numerical stability. Thus, for each ∆x j , there is a corresponding stable time step and we denote this by ∆t j . This in turn implies that the number of iterations required to estimate w j will strongly depend on j.
Once the simulations are completed and the values of w j are known, we form the three norms given by The corresponding grid convergence indices (GCIs) G j are then given by where p is the order of accuracy of the finite difference scheme. In our case, since the scheme is first-order accurate in space, we set p = 1. According to the theory presented in [66], if G j < G j+1 for j = 1 and j = 2, then the finite difference scheme is grid convergent. The finite difference simulation was performed using the Matlab ® script with the parameters presented in Table 1.
The results of the simulations are presented in Table 2 below. The norms are then given by Finally, computing the GCI values, we have: These are good results. We see that G 1 > G 2 > G 3 , which implies that the finite difference scheme fulfils the criteria for grid convergence given by Richardson [66]. We are concerned that while the values of w j demonstrate grid convergence, this analysis does not measure the number of oscillations of the beam tip before t . To verify whether the number of oscillations is stable between cases, a second grid convergence study was performed with β set to zero throughout, and the results are presented in Table 3. Table 3 includes the extra column indicating the number of oscillations performed by the vibrating beam during the simulation time. The spread of values of w j was explained by the fact that the measurement was taken at different phases of oscillation. The amplitudes of oscillation were very closely matched.   This grid convergence study was performed with the algorithms implemented in Matlab ® . On a PC with 8 Intel Core i7-7800X @3.50GHz processors, one iteration of the fourth case of the grid convergence study (101 nodes) takes approximately five microseconds to execute, as measured by the Matlab ® tic and toc functions. This way, to compute one second of beam dynamics for this case takes approximately 10 s of computation time.
In the next sub-section, we describe activities undertaken to compare the predictions using the finite difference approach with the analytical solution.

Comparison of the Simulation Results with an Analytical Solution
Here, we compare the finite difference algorithm with the analytical solution of the homogeneous equation derived in Appendix A. We compare the time histories of the tip deflection of a beam without damping predicted by both approaches for the specific boundary and initial conditions of a cantilevered beam with initial deflection equal to the first mode shape and without an external load. The geometry and structural properties used here are the same as for the beam in the grid convergence study presented above.
For the purposes of validating the finite difference algorithm against the homogeneous Euler-Bernoulli equation, we start with Equation (A28) and set the initial geometry of the beam to: with H 1 = 1. As such, we dramatically simplify the form of the function w(x, 0) without compromising the method of validation of the finite difference method. Then, we define the exact analytic function w a as w a (x, t) = (cos α 1 x − cosh α 1 x) For a beam of length L = 10 m we, find that the first positive value of α (denoted α 1 ) satisfies Equation (A16) to be 0.1875.
With E = 70.0 × 10 9 (aluminium), I = 6.67 × 10 −5 and M = 80 (beam mass evenly distributed over all the nodes), it follows that:  Figure 8 presents the deflection of the 10-m cantilevered beam in free vibration using Equation (36) as a function of length along the beam x and time t. Figure 9 presents a comparison of the tip excursions predicted for this cantilevered beam in free vibration using Equation (36) and using the finite difference algorithm presented above. It is clear that the natural frequency of vibration predicted using the finite difference approach is higher than the analytical solution when fewer nodes are used, and as the number of nodes increases, the predicted frequency decreases towards the analytical solution.

Application of the Finite Difference Scheme to a Realistic Case of a Commercial Aircraft Wing in a Gust Encounter
As part of the certification procedure for any new design of commercial aircraft, the design authority (manufacturer) is obliged to demonstrate that the aircraft can satisfactorily operate when subjected to gusts and turbulence during flight. A set of standard vertical gust velocity profiles is specified by the certification authority for use in the certification process, and one of these profiles is often referred to as the "one minus cosine" profile. A full description of the gust specification can be found in the Federal Airworthiness Administration publication [68]. The one minus cosine profile is a vector field defined as v(x, y, z) = where i, j, and k are the unit vectors in the three orthogonal axes x, y, and z, respectively, and where M and S are fixed scalar constants.
Here, we present the results of a set of simulations of the finite difference algorithm applied to the case of a wing of a large commercial passenger aircraft encountering a one minus cosine gust during steady forward flight to predict the vertical deformation of one of the wings. We assumed the aircraft is symmetric with regard to the vertical centre-line plane, and that the wings are built into a rigid fuselage; therefore, it is only necessary to analyse one wing (the right wing in this case). We estimated distributions of mass, the second moment of area, and Young's modulus from the general engineering experience, and we used an inhomogeneous source function predicted by unsteady potential flow theory for the gust encounter. Figure 10 shows the results of the simulation that starts at time t = 0 with zero initial deflection; that is, w(x, 0) = 0 for all x.
The Euler-Bernoulli Equation (1) models the transverse translation deformation of the beam. It does not explicitly incorporate any coupling with twisting deformation. In physical aero-elastic scenarios, there is strong interaction between the bending and twisting. For the case described in this section, only bending was considered. This was done to specifically study the finite difference method in isolation. This way, we were sure that the dynamical deformation of the beam observed here was only due to the algorithm presented in this paper.
An Unsteady Vortex Lattice Method (UVLM) developed by one of the authors and implemented in the software package Flexit and described in [24] was used to predict the inhomogeneous source function q(x, t). The UVLM was described in [69]. The gust velocity profile (which varies with the position in the flow field) affects the boundary conditions on the UVLM as the aircraft flies through the gust. These changes in the boundary conditions affected the resulting predicted pressure distribution over the wing. This change in pressure distribution in turn changes the Euler-Bernoulli source function q(x, t), and the finite difference method predicts a change in the deflection function w(x, t) of the wing spar. It is crucial to appreciate that this change in w(x, t) results in a change in the wing aerodynamic geometry and a new bound vortex geometry. On the next time step, the UVLM uses this new bound vortex geometry and the resulting pressure distribution prediction is based on this new geometry.
The unsteady nature of the computation ensures that there is a feedback between the fluid-dynamic forces and the wing geometry. The resulting simulation is therefore not simply a forced vibration, but rather a full time-varying fluid-structure interaction (FSI) simulation. In the present work, we used the Flexit software code to perform these calculations. More details of the FSI algorithms used in this software are described in reference [24].
For the purpose of this demonstration of the algorithm, we discretised the wing structural beam into just eight nodes. For a typical aeroelastic analysis, many more nodes would be used. The left-hand plot shows a snapshot of the deformation geometry at the end of the simulation. The right-hand plot shows the wing tip deflection and the value of the load applied along the wing as they vary throughout the simulation. The beam oscillates at the start of the simulation as the transverse load is initially applied. The amplitude of oscillation decreases, and at approximately 1 s, the oscillation has almost completely died out (and the tip deflection settles at approximately 4 m). At 2.2 s, the aircraft flies into the one minus cosine gust profile, and the potential flow theory predicts a transient increase in the transverse load. The finite difference algorithm then predicts a sudden increase in wing tip deflection peaking at approximately 8 m before again reducing to 4 m after the aircraft has flown past the gust profile. The algorithm presented in this paper was included in a computational scheme where the twist and bending algorithms were executed in parallel. Implicit coupling takes place in this scheme because at each time step, the twist of the wing changes the UVLM panel geometry. This changes the predicted aero-dynamic pressure distribution, which in turn modifies the transverse load distribution q(x, t) in Equation (5). This coupling is not the immediate topic of this paper; however, the reader can see an illustration of the effect of the twisting on bending implemented in this way, in Appendix F.

Damping Coefficient β
For the established approach to predict the deformation of a beam using modal analysis, Rayleigh damping is usually used, with a mass damping coefficient µ and a stiffness damping coefficient λ. Estimates of µ and λ were established over many years by comparison between experiment and simulation.
For the approach developed in this work, the damping term is velocity-dependent, and we did not find a quantitative correspondence between the two damping methods.
As a result, it is difficult to estimate the correct value of the damping coefficient β that should be used when performing simulations with the finite difference algorithm. Structural damping is highly sensitive to the fine details of how the beam is constructed and assembled, and it is also sensitive to temperature. However, comparing the magnitudes of the damping force with typical external applied loads, it is clear that the effects of the damping forces are relatively small. In simulations for the design of an aircraft wing, the authors noted node deformation velocities of the order of 5 m per second. For a typical value of β of 10, this results in a damping force of 50 Newtons per node, while the external load applied to the node is approximately 1000 Newtons. Thus, the damping force accounts for approximately 5 percent of the response of the beam. If the error in the estimate of β is 50 percent, for instance, then the error in the simulation that results will only be 2.5 percent.

Comparison of Static Deflection Experiment with Finite Difference Simulation
As part of a general research project being performed by the Dynamics, Control and Simulation group at Cranfield University, experiments were conducted to determine the steady-state equilibrium tip deflection of a scale model of an aircraft wing spar clamped at the wing root when different transverse loads are applied at the tip. Here, we present the results of the work we performed to compare the static deflections measured during the experiments with simulations performed using the finite difference algorithms presented in this paper. Note that here we only compare the steady-state tip deflections once any transient oscillations have died out. More examples of the comparisons of beam bending theory with experimental data can be found in [70]. Figure 11 presents the geometry and structural properties of the wing spar model. The scale model under investigation is very simple, and consists of a single flat 6082-T6 specification aluminium plate, nominally measuring 2 mm in thickness with a rectangular cross-section. Accurate measurements of the sample conducted by the authors using a micrometre showed that the actual plate thickness varied from 1.7 mm to 1.9 mm along the length of the spar. Figure 12 shows the experimental set-up. Measurements of the tip deflection were recorded with a resolution of 1 mm. It is estimated that the measurement error is +/− 1 mm (2 σ).  The wing spar model beam properties are also presented in Table 4. In the experimental set-up, the wing was clamped ad built-in over the whole rectangular area on the left of the sample, and the larger irregular area was freely cantilevered. The geometry of the experiment, the relatively small tip loadings, and the small magnitudes of the deflections imply that the Euler-Bernoulli PDE can give a good approximation of the experiment.
The dynamic simulation was performed using a number of nodes ranging from 27 to 102, a damping coefficient β of 5.0, and varying time steps (always satisfying the Lax-Richtmyer criteria), and we used the planform shape and second moment of area (I y ) distributions of the actual test piece. Plate thicknesses ranging from 1.6 mm to 1.9 mm were investigated. The simulation began with zero initial deflection and the prediction was sampled after several million iterations after the oscillations completely damped out. Table 5 presents the loading cases and the corresponding deflections measured in the experiment and predicted by computer simulations using the finite difference algorithm. The errors between the experimental results and the simulations are quite large. Figure 13 demonstrates that, as well as the magnitudes of deflection being in disagreement, there is also disagreement in the rate of the change of deflection with the load between the experiment and the simulations. Before passing judgement on the accuracy of the simulation, it is important to consider factors in the set-up of the experiment that may lead to these large discrepancies. These factors may include the impact of the inaccurate manufacture of the scale model such as variations in plate thickness, the arrangement for the clamping of the built-in spar, the method and accuracy of measurement of the deflection, and the existence of the pre-bend in the spar before the tests took place. In addition, the load was applied close to the tip, and not exactly at the tip of the beam. Figure 13 presents simulation data for a wide range of beam thicknesses in order to demonstrate that the algorithm is stable and robust over a wide range of input parameters.
We note that a small variation in the thickness of the aluminium sheet results in a large change in the predicted magnitudes of deflection. This is because the thickness is increased to the fourth power when computing the second moment of area. A change in thickness of 5 percent leads to a change in predicted deflection of 16 percent. Then, we note the error between the tip deflections predicted for zero tip loads and those measured in the test. These may be partially explained if the sample had a pre-bend before the test. Of more concern, however, is the difference in the slopes of the load versus deflection curves. The slope of the simulation predictions is lower than that for the experiment, implying that the simulation overestimates the bending strength of the beam. Finally, we note that, as the number of nodes used in the simulation increases from 27 to 102 (for a plate thickness of 1.8 mm), the error in the predicted displacement in fact increases.

Conclusions and Future Work
Here, we discuss the main conclusions concerning the finite difference method presented in this paper. The finite difference method provides an alternative approach to the established mode shape analysis for the approximation of both the homogeneous and inhomogeneous versions of the Euler-Bernoulli equation. We demonstrated that the method can be applied to the analysis of aeroelastic phenomena in aircraft. The method results in excellent agreement with the analytical solution to the homogeneous Euler-Bernoulli equation. As the number of nodes increases, the finite difference algorithm appears to rapidly converge to the analytical solution. By applying the method to an idealised commercial aircraft wing subjected to external loading, we showed that the method can be applied in the general case where E, I, µ, and ∆x arbitrarily vary with x, and where the external load q arbitrarily varies with x and t. Indeed, the finite difference algorithm is stable and well behaved when performed for this realistic physical case. The equilibrium deformation predicted by the finite difference algorithm appears to give poor agreement with a particular set of experimental results. To form an accurate picture, comparison with other experimental results is recommended. We extended the application of the Lax-Richtmyer stability criteria to the fourth-order Euler-Bernoulli partial differential equation using a finite difference approximation. The stability of the finite difference algorithm is highly sensitive to the choice of forward, backward, and central finite difference approximations to the various derivatives. This paper considered the case of a simple beam geometry with one end cantilevered. By choosing different boundary conditions, it is possible to consider beams with different anchorings. The boundary conditions at the ends of the beam can easily be controlled by setting the values of w at the ghost nodes.
In this work, we considered the deflection function w : R 2 → R that describes the scalar deflection in the z direction as a function of distance along the beam x and time t. We can extend the range of cases that can be analysed to beams subject to external load distributions and undergoing deflections in any direction in the yz plane (see Figure 1). These situations can be analysed by resolving the loads in the orthogonal y and z directions, and performing independent finite difference simulations in each of the resolved directions. The scheme described in this paper assumes that all the beam nodes are colinear. Real beams to be analysed in engineering applications are often not straight, and approximating them by a straight line is not appropriate. By expressing the Euler-Bernoulli equation in the generalised curvilinear coordinates, we can admit the analysis of beams with noncolinear nodes. The dynamical structural behaviour of a given beam undergoing angular rotation is different to that of the non-rotating beam. The effect of the angular velocity of a real beam is to increase its effective stiffness and increase the natural frequencies of vibration. By ignoring the effect of angular velocity, the predicted deformations will be too large and the oscillatory frequencies will be too low. To use the finite difference method for the analysis of rotating beams such as helicopter rotor blades, an extra term must be included in the PDE to account for centrifugal forces. See Appendix C for a description of the extra term required to take account of centrifugal forces. The effects of rotation of the cross-sections in the xz plane can become significant when larger deformations are experienced or when the shear modulus of the beam material is low. Developing a finite difference approach for the Timoshenko coupled equation will significantly extend the body of knowledge in two ways; firstly, it will provide an example of application of this method to coupled PDEs with two unknown functions of x and t (w and ϕ); and secondly, it will provide a second method of predicting beam deflections and a comparison of the two methods will contribute to their validation. A description of the Timoshenko-coupled PDEs can be found in Appendix D. Finally, the choice of biases for the finite difference approximations should be further investigated. Acknowledgments: The authors are grateful to Douglas Muirden, Geert Kapsenberg, and Michael Strauss for their informal discussions that stimulated the ideas presented here. We would also like to acknowledge the IT support at Cranfield University, UK. We would also like to acknowledge the constructive comments of the anonymous reviewers of the MDPI Aerospace journal.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The

Appendix A. An Analytical Solution of the Homogeneous Euler-Bernoulli Equation for Validation Purposes
As part of the validation of our numerical work, we will compare our approximate numerical results against an exact analytical solution, and we derive the exact solution in this section.
For the general case, we might consider a rigidly cantilevered beam at one end (and otherwise free to deform) and then a time and span-wise varying transverse loading distribution function q(x, t) is applied. However, it is difficult to derive any analytical results when q(x, t) = 0 and β = 0 in Equation (2). Therefore, purely for validation purposes, we will consider here the special case of the homogeneous Equation (q(x, t) = 0) and β = 0 with the rigidly cantilevered beam at one end. In this case, an analytical solution can be derived as outlined below.
While the homogeneous case is relatively simple, it does exercise much of the behaviour of the numerical method, and so provides a valuable means of improving confidence in the validity of the current work.
With these considerations in mind, we derived an analytical solution with boundary conditions where the beam is cantilevered at one end and the beam is free to oscillate without damping. The inhomogeneous case, in which the transverse load distribution is significant and varying in time and space will be subsequently discussed along with a comparison with experimental data.
Let us consider the homogeneous version of the Euler-Bernoulli equation with the beam property functions E, I, and µ independent of x and with the external force applied to the beam q(x, t) and the damping coefficient β both set to zero. In this case, we can write: which is often referred to as the free vibration equation. By applying the technique of separation of variables, an analytical solution to the homogeneous Equation (A1) can be found for a beam of length L of the form: where the corresponding boundary conditions are: where we adopt the notation f (p) to denote the p th derivative of the function f with respect to its independent parameter. Substituting Equation (A2) into Equation (A1), we obtain: and re-arranging Equation (A4) gives: where the left side of Equation (A5) is only dependent on x and the right side is only dependent on t. Since x and t are independent variables, both sides must be equal to the same constant, e.g., λ. Thus, we have: and: In this way, we expressed the partial differential Equation (A1) as two ordinary differential equations and a constant λ.
Beginning with the first ordinary differential Equation (A6) in x, we impose the constraint that we will only consider real-valued functions X. We justify this simplification by appealing to the fact that the deflection itself is real. We can also exclude the trivial and physically un-meaningful solutions when λ = 0 and λ < 0, respectively. This leaves us with the case λ > 0, and we consider the trial solution: where λ = α 4 . After differentiating the trial solution (A8) four times, we obtain: and: We can see that X (4) (x) = α 4 X(x) = λX(x), so the trial solution given in Equation (A8) is indeed a solution to the ordinary differential Equation (A6). It remains to determine the coefficients A, B, C, and D, and the new constant α. To do this, we use the boundary conditions (A3) of the partial differential Equation (A1). By applying boundary conditions (A3) to our trial solution (A8) and imposing the condition that α = 0, we obtain: thus: and: By combining Equations (A14) and (A15), we obtain the condition on α given by 1 + cos(αL) cosh(αL) = 0.
There are countably infinite real (and therefore ordered) solutions to Equation (A16) for α and we denote them by α n , n ∈ N. The first few solutions can then be numerically approximated. For each n, the function: X n (x) = A n cos(α n x) + B n sin(α n x) + C n cosh(α n x) + D n sinh(α n x), is a solution to the ordinary differential Equation (A6). Substituting for B n , C n , and D n using Equations (A13) and (A15), we obtain: +A n sin(α n x) − sinh(α n x) cos(α n x) + cosh(α n x) [sin(α n x) − sinh(α n x)] , and these functions X n are often referred to as mode shapes. Since the linear sum of solutions to an ordinary differential equation is also a solution, it follows that, in general, any linear superposition of the mode shape functions is also a solution. Figure A1 presents the first five modes for A n = 1, together with the total deformation shape given by 5 ∑ n=1 X n (x). Note that, here, A n is simply a scaling constant. Figure A1. First five mode shapes for a freely vibrating cantilevered beam.
We then turn our attention to the function T that must satisfy Equation (A7). Since λ = α 4 , it follows that there is a value of λ corresponding to each n, and we denote this value by λ n , and we write: Thus, Equation (A7) becomes: µ EI T (2) n + α 4 n T n = 0, (A20) thus: T Assuming that the solution must be real-valued, for each n, we consider the trial solution: where F n and G n are real-valued constants. By taking derivatives, we have: and: as required. Now, we recall that (∂w/∂t)(x, 0) = 0 for all x ∈ (0, L), so X n (x)T (1) n (0) = 0, and this forces T (1) n (0) = 0. Hence, G n = 0, and so for each n ∈ N, we have: Bringing together the expressions for X n and T n , the function w n (x, t) = X n (x)T n (t) can be written explicitly as where: and where H n = A n F n . For each n ∈ N, the function w n (x, t) satisfies the homogeneous Euler-Bernoulli partial differential equation at (A1) and the boundary conditions w(0, t) = 0 and (∂w/∂x)(0, t) = 0 ∀t > 0, and the initial condition (∂w/∂t)( For any given n, the function w n (x, 0) may or may not be equal to the initial geometry of the beam f (x). In fact: for all x ∈ (0, L) only if f (x) is a constant multiple of: To satisfy a given initial condition f (x) such that w(x, 0) = f (x), in general, it may be possible to expand f (x) as a linear sum of the mode shapes on the interval (0, L). We can express w(x, 0) as where the constants H n are chosen accordingly. Thus, the final form of the analytical solution to the homogeneous Euler-Bernoulli Equation (A1) can be written as where H n and α n for each n ∈ N are real constants. In Section 3.2, the values of w(x, t) given by Equation (A31) will be compared with approximations to w predicted by the numerical method to demonstrate the accuracy.
Thus, the expression: is a finite difference approximation for f (p) (x), accurate to q th -order in ∆x.
The coefficients a m remain to be found, and we subsequently address this challenge. We note that the conditions on the coefficients a m are If l ≤ 0 ≤ r, then by choosing r − l = p + q − 1, we see that the first matrix in Equation (A44) is square (since the column corresponding to m = 0 was removed), and Equation (A44) is a system of p + q − 1 linearly independent equations in the p + q − 1 unknowns a m .
Then, we consider the value of l (and hence r). If p + q − 1 is even, we can choose: and: r = 1 2 (p + q − 1).
This choice of l and r corresponds to a central finite difference. Then, to obtain the finite difference approximation for f (p) (x), it only remains to solve the system of equations in Equation (A44) using linear algebraic methods.
Since m ranges from −1/2(p + q − 1) to 1/2(p + q − 1), the substitution of the values of a m into Equation (A42) gives the central finite difference approximation for f (p) (x).
The estimate can be left, biased, or right by decreasing or increasing the value of l, respectively, (while satisfying the inequality l ≤ 0 ≤ r).

Appendix C. Accounting for Centrifugal Force in the Euler-Bernoulli Equation
The Euler-Bernoulli equation with the inclusion of the effect of centrifugal force is given in [23] as where T(x i ) is given by where L is the radius of the tip of the beam, R is the radius of the root, Ω is the angular velocity of the beam, and u is a dummy variable of integration. By including the finite difference form of the third term on the right-hand side of Equation (A47) in the finite difference algorithm, we will be able to account for angular velocity, and hence, we will extend the range of cases for which the finite difference method can be used. Work has already been completed to derive the finite difference form of the centrifugal term, and initial computer experiments to analyse the behaviour of a helicopter rotor blade show a good comparison with experimental data from full-scale tests-as can be seen in [24].
Bernoulli equation discussed in this paper. The script implements boundary conditions for a cantilevered beam with non-constant material properties and mass distribution that represents the scale model of the aircraft wing used in the experiments described in Section 3.4.2 of the paper. A "one minus cosine" time-dependent external load distribution is considered. In the interest of simplicity, just eight equally spaced beam nodes with lumped masses are considered. This simulation discretizes the smoothly varying external load distribution into 317 time slices, each one lasting for 1 −2 seconds in duration. The time step used in the finite difference calculation is computed using the Lax-Richtmyer criteria.

3.
massR.mat: This file contains the masses (in Kg) of each of the eight nodes used in the simulation.

4.
IyRnew.mat: This file contains the distribution of second moment of area perpendicular to the bending axis of the wing beam cross-section. The values are sampled at the locations of the eight nodes along the length of the beam.

5.
IyRnew.mat: This file contains the values of the externally applied loads (the "one minus cosine" distribution) at each node at each of the 317 time slices. 6.
calcDt.m: This file contains a Matlab ® function that is called by EBFiniteDifference.m and computes the maximum time step. For some runs of the simulation, the Lax-Richtmyer time steps are over-ridden and smaller time step values are used.

Appendix F. Coupling between Bending and Twisting
For illustrative purposes, here we present a comparison of wing tip deflection for a generic design showing simulations with and without twist in Figure A2. These simulations were performed with a value of S = 40 and M = 12 in Equation (37), and the wing encountering the gust centred at frame 200 of the simulation (corresponding to x 0 = 200). The coupling of bending with twisting in this way will be the topic of future research by the authors.