Noise-robust modal parameter identification and damage assessment for aero-structures

Abstract

Purpose

Ground vibration testing is critical for aircraft design and certification. Fast relaxed vector fitting (FRVF) and Loewner framework (LF), recently extended to modal parameter extraction in mechanical systems to address the computational challenges of time and frequency domain techniques, are applied for damage detection on aeronautically relevant structures.

Design/methodology/approach

FRVF and LF are applied to numerical datasets to assess noise robustness and performance for damage detection. Computational efficiency is also evaluated. In addition, they are applied to a novel damage detection benchmark of a high aspect ratio wing, comparing their performance with the state-of-the-art method N4SID.

Findings

FRVF and LF detect structural changes effectively; LF exhibits better noise robustness, while FRVF is more computationally efficient.

Practical implications

LF is recommended for noisy measurements.

Originality/value

To the best of the authors’ knowledge, this is the first study in which the LF and FRVF are applied for the extraction of the modal parameters in aeronautically relevant structures. In addition, a novel damage detection benchmark of a high-aspect-ratio wing is introduced.

Citation

Dessena, G., Civera, M., Pontillo, A., Ignatyev, D.I., Whidborne, J.F. and Zanotti Fragonara, L. (2024), "Noise-robust modal parameter identification and damage assessment for aero-structures", Aircraft Engineering and Aerospace Technology, Vol. 96 No. 11, pp. 27-36. https://doi.org/10.1108/AEAT-06-2024-0178

Nomenclature

Definitions Acronyms and abbreviations

BeaRDS

= Beam Reduction Dynamic Scaling;

EMA

= Experimental modal analysis;

FRF

= Frequency response function;

FRVF

= Fast relaxed vector fitting;

GVT

= Ground vibration testing;

HAR

= High aspect ratio;

LF

= Loewner framework;

LTI

= Linear time-invariant;

MAC

= Modal assurance criterion;

MOR

= Model order reduction;

N4SID

= Numerical algorithms for (4) subspace state space system identification;

OMA

= Operational modal analysis;

SI

= System identification;

SIMO

= Single-input multi-output;

SHM

= Structural health monitoring;

UAS

= Unmanned aerial systems;

VF

= Vector fitting;

XB-2

= eXperimental BeaRDS-2;

ζ_n

= Damping ratios;

μ

= Average or arithmetic mean;

ρ

= Density;

φ_n

= Mode shapes;

σ

= Standard deviation; and

ω_n

= Natural frequencies.

Introduction

System identification (SI) is a well-established field (Ljung et al., 2020) with broad applications in various engineering domains. In structural dynamics, SI is primarily used to model systems by determining their modal parameters from experimental or operational data (Mugnaini et al., 2022), a process known as modal analysis. Modal analysis includes two subdomains: Experimental modal analysis (EMA) (Dessena et al., 2022a) and operational modal analysis (OMA) (Sibille et al., 2023). OMA, particularly prevalent in civil engineering (Zanotti Fragonara et al., 2017), involves extracting modal parameters from output-only data and is also used in aeronautics for taxi vibration testing (Al-Bess and Khouli, 2024). Nevertheless, this work focuses solely on EMA approaches with the scope of identifying modal parameters, such as natural frequencies (ω_n), damping ratios (ζ_n) and mode shapes (φ_n). These parameters find applications in vibration-based damage detection (Civera et al., 2022) and model updating (Dessena et al., 2024a, 2024c). This work focuses on the former.

In terms of vibration-based damage detection, damage, defined as a change affecting the operational capability of the system (Farrar et al., 2001), is detected through structural health monitoring (SHM), using a statistical pattern recognition strategy involving operational assessment, data acquisition, feature selection and damage evaluation. Modal parameters obtained from EMA, like ω_n and φ_n, are crucial for vibration-based SHM (Rytter, 1993), with ω_n primarily assessing damage severity (Fan and Qiao, 2011) and φ_n supporting damage localisation. On the other hand, ζ_n is unsuitable as a damage indicator due to its dependence on non-structural factors (Civera et al., 2021a). For a general SHM overview, refer to Sohn et al. (2004); for vibration-based SHM using modal data, see Rytter (1993).

In aeronautics, EMA plays a critical role in ground vibration testing (GVT), an essential procedure in civil aircraft design and certification of large aircraft (Lubrina et al., 2014), as well as in unmanned aerial systems (UAS) design (Olejnik et al., 2022). GVT can be conducted on entire vehicles (Prananta et al., 2016) or specific components like wings or helicopter blades (Weber et al., 2021), facilitating the validation of finite element models used in design processes (Promio et al., 2018).

Current methods for extracting modal parameters face computational challenges. Issues arise from the ill-conditioning of fitting processes in the frequency domain (Lefteriu and Antoulas, 2009) and difficulties in handling large datasets in the time domain (Dessena et al., 2023a). To address these challenges, the single-input multi-output (SIMO) techniques Loewner framework (LF) (Dessena et al., 2023a) and fast relaxed vector fitting (FRVF) (Civera et al., 2021a) were recently extended to the extraction of modal parameters from frequency domain data in mechanical systems. Hence, the scope of this work is threefold:

• Compare the LF and FRVF for noise robustness and damage detection application.

• Pioneeringly apply the LF and FRVF to vibration-based damage detection in aeronautical structures.

• Introduce a novel experimental benchmark case study of a high aspect ratio (HAR) wing with simulated damage.

Firstly, the robustness to noise is assessed on an original numerical dataset for increasing noise levels. The same dataset also features different damage and loading scenarios for the preliminary assessment of the vibration-based damage detection capability of LF and FRVF on aeronautically inspired synthetic structures. Then, the original experimental dataset of a HAR wing with simulated damage is introduced. For the numerical model, the identified parameters are compared to the numerical results and for the newly introduced experimental dataset, numerical algorithms for (4) subspace state space SI (N4SID) (Van Overschee and De Moor, 1994) are used as a benchmark. The performance of the methods is assessed in terms of precision of the identification and computational performance (time-to-identification in s).

The experimental data for the GVT of the HAR wing are made available to interested readers in an openly available Zenodo repository entry (Dessena, 2024).

Methods

A thorough review of the SI discipline, or its use for modal analysis, is not within the scope of this work and the interested reader is referred to the book (Ljung, 1987) for SI, and to the classical work (Maia, 1988) for the SI role in modal analysis. Nevertheless, an introduction to the two recently introduced SI methods used within this work, LF and FRVF, is given.

Fast and relaxed vector fitting

The FRVF, introduced in Deschrijver et al. (2008), refines the vector fitting (VF) algorithm (Gustavsen and Semlyen, 1999) for modelling large multiport electrical circuits. FRVF enhances VF by introducing a relaxed constraint during pole identification (Grivet-Talocia and Gustavsen, 2016; Gustavsen, 2006) and using QR decomposition (such that a matrix A¯ = QR of an orthonormal matrix Q and an upper triangular matrix R) for efficient matrix computation. In the form considered here, the FRVF procedure was tested and applied for the first time for the SI of simple mechanical systems in Civera et al. (2021b) and to large civil structures and infrastructures in Civera et al. (2021a).

FRVF is an iterative method. Following a linear least-squares approach, it relocates initial poles to minimise the divergence between estimated transfer functions and experimentally recorded data. That is to say, the complex-valued data f(s) ∈ C^p×k (for k frequency samples, defined over s = jω and for p frequency response functions [FRFs]) are approximated (for a generic single input-multiple outputs configuration) by the rational function:

(1) f(s)≈ ∑ m=1 k c m s− a m +D+se ^

where k is the preset number of poles to be identified, a_m is the mth pole (at any given iteration), c_m is the corresponding vector of residues, D ∈ R^p×1 is the feedthrough matrix (again, reduced to a column vector for a single input scenario) and e^ ∈ R^p×1 contains the terms proportional to s (generally all zeroed and included here for completeness only). The poles are all complex conjugates, thus they come in pairs as a_m = −α + jβ, (a_m + 1) = −α − jβ, where the real part is strictly negative (i.e. they all lie in the left half of the complex plane) to ensure the stability of the dynamic system. For the first set of poles, β is linearly spaced over the frequency range of interest and α = β/100. These represent the departing points for the iterative, deterministic optimisation procedure.

The components of the rational model described in equation (1) can then be seen as the parameters of the state-space model, such as:

(2) f(s)≈C (sI−A) −1 B+D+se ^

where, according to the common terminology of state-space modelling, A ∈ C^k×k is the system coefficient matrix, B ∈ R^k×1 is the input matrix (which, for this application, is reduced to a column vector of ones) and C ∈ C^p×k is the output matrix. The VF identifies the poles of f(s) by solving (in the least-square sense) the following linear problem:

(3) σ(s)f(s)=p(s)

where:

(4a) σ(s)= ∑ m=1 k c ˜ (s) s− q m +1

(4b) p(s)= ∑ m=1 k c ˜ (s) s− q m +D+sE

indicate the mth tentative pole. In the fast and relaxed variant of the VF variant, the definition of σ(s) is slightly changed; specifically, it is relaxed as:

(5) σ(s)= ∑ m=1 k c ˜ m s− q m +D ˜

where D¯ is real but not necessarily unity (even if it is expected to approach 1 as the procedure converges throughout the subsequent iterations). Thus, the pole relocation procedure aims at achieving the terms p_m that most closely resemble their (unknown) counterparts (a_m), such that the root-mean-square error of approximation for f(s), i.e. the error of the fitting, is minimised.

More details concerning the implementation and the several technical aspects of FRVF can be found in (Grivet-Talocia and Gustavsen, 2016).

Loewner framework

The LF was first introduced as a SIMO SI technique for the extraction of modal parameters for vibration-based SHM of mechanical systems in the frequency domain in Dessena et al. (2023a).

The LF was notably proposed in Gosea and Antoulas (2015) as a model order reduction (MOR) technique for large dynamical systems, but its origins can be traced back to the 1930s, when Charles Loewner introduced the namesake interpolation matrix (L) (Löwner, 1934).

Antoulas et al. (2017; 2007) have developed the LF for the MOR of dynamical systems by considering tangential interpolation, or rational interpolation along tangential directions (Kramer and Gugercin, 2016). Later, they applied it to the SI of electronic systems in Lefteriu and Antoulas (2009) to relax the severely ill-conditioning of current fitting processes (Lefteriu and Antoulas, 2010).

Given a linear time-invariant dynamical system ∑ with k internal variables in time-continuous descriptor-form representation, m inputs and p outputs:

(6) ∑:E d dt x(t)=Ax(t)+Bu(t);y(t)=Cx(t)+Du(t)

where x(t) ∈ R^k×1 is the internal variable, u(t) ∈ R^m×1 is the function input and y(t) ∈ R^p×1 is the output. The constant matrices for k internal variables are:

(7) E,A∈ ℝ k×k ,B∈ ℝ k×m C∈ ℝ p×k D∈ ℝ p×m

When, for a given finite value λ, the matrix A- λE is non-singular, such that λ ∈ C, a Laplace transfer function, H(s), of Σ can be defined in the form of a p by m rational matrix function:

(8) H(s)=C (sE−A) −1 B+D

The LF, via tangential interpolation, actively fits the FRF data to H(s). The aim does not differ from more established techniques, such as the rational fraction polynomial method, but the definition of the tangential interpolation directions (random in practice) allows for a more streamlined computation. A theoretical exposition of the LF is found in Mayo and Antoulas (2007), while the recent contributions to the implementation of the LF are outlined in full detail in Dessena et al. (2023a), a comparison of its computational performance is found in Dessena et al. (2023b) and a MATLAB tutorial is given in Dessena (2023a, 2023b).

Numerical case study

For comparing the robustness to noise of the LF and FRVF, a numerical case study of a wing spar of a small (maximum take-off mass <7 kg) fixed-wing UAS is introduced. The spar is discretised as a four-element two-dimensional Euler–Bernoulli cantilever beam in MATLAB by using standard (4 × 4) mass and stiffness matrices [see equations (A1a) and (b) in the Appendix], where the degrees of freedom are rotations (θ^xy) and vertical displacements (x_n), as shown in Figure 1(a). The spar square box cross-section is shown in Figure 1(b).

The beam is made from 6061-T6 aluminium with an elastic modulus (E) estimated at 70 GPa and a density (ρ) of 2,700 kgm⁻³. The beam section dimensions are defined in Figure 1(b), while its length (L) is 1.1 m.

For comparing the SHM capabilities of LF and FRVF, five different scenarios exist for the spar:

1. baseline condition;

2. 5% stiffness reduction in the third element;

3. 10% stiffness reduction in the third element;

4. 30% stiffness reduction in the third element; and

5. Pylon and engine at the spar’s midpoint, discretised as a 0.3 kg lumped mass.

These were selected to have four progressive damaged scenarios and a change in mass configuration, which could mimic an under-wing payload or a pylon-engine assembly. The ζ_n is set to 3% for all cases and modes.

Investigation of noise effects

The LF and FRVF robustness to noise is tested by corrupting both input and output channels with additive white Gaussian noise at 0%, 0.1%, 0.2%, 0.5%, 1%, 1.5%, 2%, 3%, 4% and 5% of the signal standard deviation (σ), totalling ten independent cases per scenario. Because LF is a numerical method primed by random starting points, the tangential directions which are random in practice (Quero et al., 2019), a numerical study over the ten noise cases for the five scenarios is carried out by running the identification 100 times (100 realisations) at the minimum order. This means that modes appearing in at least 90% of the realisations are considered. On the other hand, FRVF is a fully deterministic method; therefore, such a study is not needed. Nevertheless, FRVF is an iterative process; thus, for the scope of this work, the number of iterations was set to 5 at its minimum order. The minimum order of k for detecting eight modes for both LF and FRVF is 16 (e.g. number of systems poles).

Figures 2(a) and (b), compare the identified ω_n from LF and FRVF and the numerical values. The modes are considered identified if less than a 10% deviation from the expected numerical result is found. For brevity, only the results of scenarios # 1 and 5 are graphically presented as the behaviour in scenarios # 2–4 is similar to # 1.

In Figure 2(a), the LF demonstrates superior capability in identifying a wider range of modes than FRVF. However, LF struggles with higher modes as noise levels increase, unlike FRVF, which consistently identifies the first mode but is less reliable with modes 2–5. LF generally identifies correctly the ω_1–5 across noise levels, with perfect identification in noiseless and 0.1% noise scenarios. Conversely, FRVF accurately identifies all modes in the first three noise cases except for the fifth scenario, where the first mode is correctly identified only in the noiseless case due to initial pole convergence challenges at lower frequencies when artificially added noise is present.

Concerning the ζ_n (shown as damping in %) of the stable modes, Figures 2(c)–(e) (LF) are more coherent with the expected values, when compared to Figures 2(d)–(f) (FRVF). LF consistently approaches the expected value across noise scenarios, while FRVF shows full consistency only in no and 0.1% noise cases, with partial coherence in the 0.2% case. The comparison of mode shapes φ_n via the modal assurance criterion (MAC) (Allemang and Brown, 1982) shows LF [Figures 2(g) and 2(i)] and FRVF [Figures 2(h) and 2(j)] results largely coherent with expected values, with LF showing more consistency overall.

In summary, the LF struggles with higher modes. However, the lower modes are identified inconsistently. Compared to the FRVF, the LF tends to produce more coherent results for damping ratios and mode shapes across noise levels.

Investigation of damage effects

This section compares LF and FRVF for damage detection using noiseless data, focusing on ω_n and φ_n, while ζ_n remains constant. Figure 3(a) illustrates the percentage difference between ω_1-3 identified by LF and FRVF for damaged cases relative to the baseline, confirming the sensitivity of both methods to damage detection, even for small damage scenarios. In addition, Figure 3(a) demonstrates close agreement between identified and numerical values. LF and FRVF effectively detect and localise small changes, with the largest deviations (in terms of modal displacement) occurring at the damaged elements. φ₁ deviations increase with damage and mass addition, especially noticeable between nodes 3 and 4 [Figures 3(b)–(c)], demonstrating the suitability for both LF and FRVF SHM.

With regards to the time to identification, the performance of the two methods, for the minimum order k = 16, is reported. The mean (μ) and σ for the time to identification, are, respectively, 4.170 s (LF) and 0.375 s (FRVF) and 0.381 s (LF) and 0.264 s (FRVF). The statistical measures were computed over 50 identification realisations. The results are obtained from a Windows 11 desktop machine with an Intel® Xeon® Processor E5-1650 at 3.20 GHz and 32 GB of RAM. From the results, FRVF is both the least computationally demanding method and the most stable in time, because its σ is smaller than that from LF.

In conclusion, LF and FRVF succeed in identifying modal parameters under no or low noise conditions, with LF demonstrating greater resilience to higher noise levels. LF accurately predicts low and medium-frequency modes, while FRVF does better with high-frequency modes. LF robustness to noise enhances its efficiency in real-world applications, particularly for detecting structural changes like damage or added mass. Although FRVF outperforms LF in computational efficiency, the former identifies more stable modes across scenarios, indicating better identification capabilities.

Ground vibration testing of a flexible wing

After the successful implementation of LF and FRVF on a synthetic benchmark, a new experimental case study, based on the eXperimental BeaRDS-2 (XB-2) HAR wing (Figure 4) developed within the Beam Reduction Dynamic Scaling (BeaRDS) project at Cranfield University (Hayes et al., 2019; Pontillo, 2020; Pontillo et al., 2018; Yusuf et al., 2019), is introduced.

The XB-2 wing was conceived as a dynamically scaled example of a civil jet airliner wing to be tested in the university wind tunnel. The wing is made of three components: the spar, the stiffening tube and the skin.

The aerodynamic surface of the wing, outlined by a NACA 23015 aerofoil, has a span of 1.5 m (1.385 m from the origin in Figure 4), with a mean aerodynamic chord of 0.172 m, a taper ratio of 0.35, a leading-edge sweep of 1.49° and a mass of 3.024 kg. The wing has neutral twist and dihedral angles. The skin is responsible for transferring the aerodynamic loads to the underlying structure and is made of two 3D-printed plastics: Digital ABS and a rubber-like compound, Agilus 30. Combining a rigid and a rubber-like material allows the skin to be flexible and preserves structural integrity.

The wing torque box consists of the spar and tube assembly. The spar was machined from two 6082-T6 aluminium blocks which were welded together and secured with four bolted L-profile plates. The main spar features a Saint George cross-shaped cross-section and a variable taper along its span, while the tube is a stainless steel tube with a 10 mm diameter and 1 mm thickness.

For the scope of this work, four scenarios are considered. The specimen described above is considered the baseline and three loaded scenarios are introduced, as described below:

• Baseline: Total mass 3.024 kg;

• Added masses: 75 g at 1,010 mm, 12 g at 1,050 mm and 61 g at 1,365 mm – 3.172 kg;

• Added masses: 88 g at 1,010 mm, 51 g at 1,050 mm, 83 g at 1,205 mm and 61 g at 1,365 mm – 3.307 kg; and

• Added masses: same as Scenario 2 plus 181 g at 570 mm and 170 g at 665 mm – 3.658 kg.

The scenarios represent progressive addition of mass to simulate increasing damaged scenarios on the wing. The location and magnitude of the added mass are limited by the existing configuration of the wing, such as ballast mass fixing points and removable panels location.

The three loaded scenarios are introduced to assess LF and FRVF sensitivity to changes in aeronautical structures. In addition, the LF and FRVF results are compared with modal parameters extracted with the well-established method N4SID. From a previous modal survey (Dessena et al., 2022b) involving the baseline wing, the first three dominant modes in the vertical direction were identified between 3 and 20 Hz. Hence, for this work, the frequency band for the linear sine sweep excitation was set between 2 and 25 Hz and spanned across 20 min (single sweep). The data is available in the data availability statement. A thorough description of the hardware setup is found in Dessena et al. (2022b).

Stabilisation diagrams help identify modal parameters from the experimental data, and Table 1 presents the results for the stable modes. LF and N4SID generally had the most stable diagrams, but even so, this did not undermine the precision of FRVF.

In Table 1, the stable ω_n and ζ_n identified with N4SID, LF and FRVF are presented. The identifications with LF and FRVF are mostly consistent with the N4SID benchmark results. In fact, the relative difference between the ω_n of LF and FRVF and those from N4SID never exceeds 1%, while the ζ_n is mostly under 10% and on rare occasions between 10% and 15%. However, this is expected due to the intrinsic nature of damping (Civera et al., 2021b).

In general terms for φ_n, the diagonal values of the MAC matrices – the off-diagonal terms have a negligible value – between the identified values from LF, FRVF and N4SID, for all cases, apart from φ₄ in the baseline scenario for FRVF, exceed 0.95 (MAC = 0.89 for φ₄ of the baseline scenario), which means an almost perfect correlation between the modes and presented methods. This is shown in Figure 5.

Figure 6 shows the φ₁ identified by the three methods in the baseline and scenario # 4. A deviation from the baseline φ₁ is displayed for all methods. In addition, a change is shown for all relevant positions (the three spanwise points corresponding to the accelerometer locations), accounting for the simulated damage localisation. Only the baseline and # 4 scenarios are shown for conciseness and clarity, but similar results are available for the other scenarios.

Conclusion

This study compares the LF and FRVF for modal parameter extraction from vibration data. Key findings include:

• FRVF is more affected by measurement noise but faster than LF by an order of magnitude.

• FRVF detects higher-order modes better, while LF identifies more modes accurately with precise damping estimation.

• LF provides more stable results with fewer spurious identifications.

• Both align with standard methods (N4SID), with FRVF slightly more accurate and LF more stable.

• Both detect damage and structural changes accurately.

• Furthermore, a novel benchmark for high-aspect-ratio wing damage assessment is introduced.

Acknowledgements

The authors wish to thank Dr Ivan Petrunin from the Centre for Autonomous and Cyber-Physical Systems at Cranfield University for providing the facilities and part of the equipment used for the tests and Prof Stefano Grivet-Talocia of the Department of Electronics and Telecommunications of the Politecnico di Torino for his valued advice.

Funding: The authors from Cranfield University disclosed receipt of the following financial support for the research, authorship and/or publication of this article: This work was supported by the Engineering and Physical Sciences Research Council (EPSRC) (grant number 2277626). The second author is financially supported by the Sustainable Mobility Center (Centro Nazionale per la Mobilità Sostenibile – MOST), Spoke 7 (Cooperative Connected and Automated Mobility and Smart Infrastructures), Work Package 4 (Resilience of Networks, SHM and Asset Management)

Author contributions: Conceptualisation, G.D. and M.C.; methodology, G.D., M.C. and A.P; software, G.D., M.C. and L.ZF.; validation, G.D., M.C., A.P., D.I.I. and J.F.W.; formal analysis, G.D. and M.C.; investigation, G.D.; resources, A.P., D.I.I., J.F.W. and L.ZF.; data curation, G.D. and M.C.; writing–original draft preparation, G.D.; writing–review and editing, G.D., M.C., A.P., D.I.I., J.F.W. and L.ZF.; visualisation, G.D., M.C. and D.I.I.; supervision, D.I.I., J.F.W. and L.ZF.; funding acquisition, L.ZF.

Conflict of interest: The author(s) declare none.

Data availability statement: Experimental data underlying this study can be accessed through the Zenodo Repository at https://doi.org/10.5281/zenodo.11635814 under the terms of under the terms of [GNU General Public License (GPLv3)].

Figures

Figure 1

4-element beam: beam section and profile

Figure 2

4-element beam: ω_n identified for scenarios # 1 and 5

Figure 3

4-element beam: change in the identified ω_1-3 and φ₁ wrt the damage scenario

Figure 4

XB-2 wing top view. y-axis coming out of the page

Figure 5

XB-2 wing: φ_1,2,4 identified by LF, FRVF and N4SID for the baseline scenario

Figure 6

XB-2 wing: φ₁ of the baseline and Scenario #4 identified by LF, FRVF and N4SID

Table 1

Natural frequencies and damping ratios identified by LF and FRVF for all scenarios

Mode	1st bending			1st coupled			2nd coupled
Scenario	N4SID	LF	FRVF	N4SID	LF	FRVF	N4SID	LF	FRVF
Natural frequency [Hz]
1	3.190	3.202	3.203	11.896	11.886	11.858	17.763	17.703	17.725
2	2.957	2.958	2.945	12.096	12.134	12.083	17.350	17.302	17.294
3	2.775	2.769	2.788	12.002	12.025	12.014	17.079	17.101	17.023
4	2.729	2.725	2.727	11.970	11.965	11.938	15.067	15.052	15.004
Damping ratio [−]
1	0.032	0.040	0.028	0.066	0.063	0.065	0.058	0.061	0.062
2	0.021	0.024	0.025	0.060	0.057	0.058	0.061	0.056	0.060
3	0.019	0.022	0.021	0.058	0.055	0.057	0.050	0.050	0.057
4	0.019	0.021	0.019	0.050	0.048	0.052	0.046	0.039	0.038

Source:

Created by authors

Appendix

The element mass and stiffness matrices for the 2D Euler Bernoulli beam considered are presented below. L is intended here as the length of a single element.

(A1a) M e = ρAL 420 [ 156 22L 54 −13L 22L 4 L 2 13L −3 L 2 54 13L 156 −22L −13L −3 L 2 −22L 4 L 2 ]

(A1b) K e =[ 12EI L 3 6EI L 2 − 12EI L 3 6EI L 2 6EI L 2 4EI L − 6EI L 2 2EI L − 12EI L 3 − 6EI L 2 12EI L 3 − 6EI L 2 6EI L 2 2EI L − 6EI L 2 4EI L ]