Quantification of Uncertainty in the Mathematical Modelling of a Multivariable Suspension Strut Using Bayesian Interval Hypothesis-Based Approach

Mathematical models of a suspension strut such as an aircraft landing gear are utilized by engineers in order to predict its dynamic response under different boundary conditions. The prediction of the dynamic response, for example the external loads, the stress and the strength as well as the maximum compression in the spring-damper component aids engineers in early decision making to ensure its structural reliability under various operational conditions. However, the prediction of the dynamic response is influenced by model uncertainty. As far as the model uncertainty is concerned, the prediction of the dynamic behavior via different mathematical models depends upon various factors such as the model's complexity in terms of the degrees of freedom, material and geometrical assumptions, their boundary conditions and the governing functional relations between the model input and output parameters. The latter can be linear or nonlinear, axiomatic or empiric, time variant or time-invariant. Hence, the uncertainty that arises in the prediction of the dynamic response of the resulting different mathematical models needs to be quantified with suitable validation metrics, especially when the system is under structural risk and failure assessment. In this contribution, the authors utilize the Bayesian interval hypothesis-based method to quantify the uncertainty in the mathematical models of the suspension strut.


Introduction
Since the mid-1960s, mathematical simulation models have been used as a tool in the field of scientific research as well as in the design and development of engineering systems [1]. Mathematical models are used to predict the non-deterministic dynamic response of engineering systems under various static and kinetic loading conditions [2]. The predicted responses are used by engineers to evaluate the system's stability, strength and vibrational behaviour under various operational conditions. However, the mathematical models are influenced by the uncertainty in the assumed model parameters, the initial and boundary conditions, the input variables as well as the governing functional relations between the model input and output variables [3]. Hence, the resulting mathematical models that predict the dynamic responses differently need to be assessed with suitable validation metrics [4], [5], [6]. A validation metric gives a quantitative measure of agreement between the predicted response and the experimentally observed response [7], [8]. The validation metric must take into account the uncertainty in the mathematical models as well as in the experimental observations of the dynamic response of the engineering system in interest [9].
The validation metrics are used for model selection when alternative mathematical models are being considered. Various validation metrics such as the classical hypothesis testing, Bayes factor, frequentist's metric and area metric are discussed in [10] and [11]. As for the probabilistic methods, validation metrics such as B factor which is based on B probabilistic approach can be utilized to judge the capability of prediction of a given mathematical model with the experimental data [12]. Model validation is classified as either univariate or multivariate. Bayesian validation metrics for univariate problems have been discussed in [13], [14]. However, for multivariate engineering systems under uncertainty different validation metrics based on the Bayesian hypothesis appproach and the area metric approach are discussed in [15], [16]. Data uncertainty occurs when the input parameters, such as state variables or the system parameters of the mathematical models are uncertain. This can be represented using probabilistic and non-probabilistic approaches, [17], [18] and [19]. The sources of uncertainty can be broadly classified as aleatoric or epistemic [10]. In this contribution, the authors consider the B interval hypothesis as a model validation metric according to [20] to quantify the uncertainty in the mathematical modeling of a suspension strut. However, and as a new approach, the authors present a way to adequately compare different models based on linear, nonlinear, axiomatic or empiric assumptions of functional relations. This work considers a suspension strut which is referred as the modular active spring-damper system (German acronym MAFDS). MAFDS is developed under the framework of German Collaborative Research Center SFB 805 at the Technische Universität Darmstadt. The concept of MAFDS is registered as a Patent DE 10 2014 106 858.A1 "Load transferring device" [21]. Fig. 1 (a) illustrates the MAFDS. The MAFDS serves as an academic demonstrator to investigate the different methods and technologies to control uncertainty in all development stages such as the conceptual phase, design and optimization phase, production and assembly phase as well as the final operational phase. The MAFDS has similar dynamic requirements of a typical suspension strut such as an aircraft landing gear. However, it is not designed to be used as a substitute to an existing aircraft landing gear. Rather, the MAFDS is an academic example to study the uncertainty in the dynamic outputs of suspension strut systems when different passive, semi-active and active modules are used to lower the uncertainty and to enhance the performance. It carries a payload on an upper truss that is connected to the lower truss via guidance links. The guidance link enables the kinematic motion between the lower and upper truss. The elastic foot acts as an impact absorption element which introduces the axial and lateral forces due to impact into the MAFDS.

Modular Active Spring Damper System MAFDS
The axial forces are transferred to the upper truss via the spring-damper component that is loaded by axial loads only due to load redirection in the upper truss of the MAFDS, Fig. 1 (b). The load redirection in the upper truss is mentioned as a key function in the patent of MAFDS "Load transferring device" [21]. An impact scenario is created by releasing the MAFDS from a specified drop height h. Uncertainty is investigated especially for MAFDS's main operating function and purpose: sustain stability, well balanced load distribution and ability to attenuate vibrations. Passive, semi-active and active technology approaches for stability, load distribution and vibration control are developed to compensate uncertainty in the main functions mentioned above. They are integrated in a modular way within the SFB 805 [26], [27], [28], [19], [22], [23], [24] and [25].
In this contribution, the global dynamic outputs of the MAFDS are the relative compression and the forces exerted by the spring damper component, and the ground impact force at the elastic foot under drop test scenario. For this, uncertainty as the varying payload and varying drop heights are assumed numerically via different mathematical models subject to experimental validation via test rig.

4
Uncertainty in Mechanical Engineering III Similarly, the elastic foot is represented by the stiffness parameter k ef . The 2DOF mathematical model has two degrees of freedom that are defined by the local translational displacements z u and z l . The 2DOF mathematical model in Fig. 2 is assumed to be under free and full homogeneous field of acceleration due to gravity g. The model is raised to the drop height h for the free fall drop test in order to excite the 2DOF-model in Fig. 2. Depending upon the drop height h, the model is excited by the initial conditions that are applied to m u and m l . Thus the initial displacements z 0 and the initial velocitiesż 0 at time t = 0 are and hence the velocityż(h) due to drop height h becomeṡ The goal is to predict the maximum dynamic response vector that is related to the displacements z u , z l and the velocitiesż u ,ż l . It consists of the maximum modeled compression z r,m,max , maximum spring-damper force F sd,m,max and the maximum ground impact force F ef,m,max .
Applied Mechanics and Materials Vol. 885 5

Four Different Approaches to Model the Stiffness and the Damping
The 2DOF mathematical model can be described with four underlying assumptions with respect to the model parameter's functional relations: linear or nonlinear, axiomatic or empiric. The authors present four different approaches to model the stiffness and the damping of the 2DOF mathematical model for the prediction of the dynamic response y vector of the 2DOF model.

a) Model with linear axiomatic stiffness k and linear axiomatic damping b
According to Fig. 2 (b), the free body diagram forces are with the linear axiomatic viscous damping force F b , the linear axiomatic spring force F k and the force F kef exerted via the elastic foot. The governing homogeneous differential equation of the 2DOF mathematical model for m u is and for m l The matrix form representation of (5) and (6) is expressed as with the global mass, damping and stiffness matrices M , B and K.

Uncertainty in Mechanical Engineering III b) Model with nonlinear empiric stiffness k(z r ) and linear axiomatic damping b
In this mathematical model, the stiffness k(z r ) of the spring is not constant and is modelled as a function that is dependent upon the value of the relative displacement z r = (z u − z l ). The empiric stiffness data points in Fig. 3 are obtained from the experimental static relative displacement under varying static loads that is conducted via the MAFDS test rig. The relative displacement z r increases nonlinearly for z r > 0.07 m. The spring force is modelled as and with the cubic stiffness parameter k 2,2 = 6.5 · 10 6 Nm −3 , quadratic stiffness parameter k 2,1 = 4.9 · 10 5 Nm −2 and linear stiffness k 2,0 = 3.3 · 10 4 Nm −1 . As seen in Fig. 3, the nonlinear regression has only a minor slope in k 2 (z r ) that is small compared to the linear regression assumption for k 1 . The matrix form representation of the homogeneous governing differential equation is expressed as

c) Model with linear axiomatic stiffness k and nonlinear empiric damping
The damping forces F b in the mathematical models a) and b) are based on linear assumptions. The damping force F b (ż r ) is obtained via experimental drop tests of MAFDS. It is seen that the empiric damping force F b (ż r ) that is obtained from the experimental test of MAFDS is related to the relative velocityż r nonlinearly. Fig. 4 (a) shows the compression time signal that is recorded experimentally by drop tests. The noise in the compression signal was filtered with a lowpass B filter. The relative velocityż r shown in Fig. 4 (b) is calculated from the measured compression z r via numerical differentiation over time. Eventually, F b (ż r ) is seen in Fig. 4 (c) as a function ofż r . The empiric damping force F b (ż r ) can be modeled as with the nonlinear damping behavior described by an adapted regression function with the cubic damping coefficient b 3 = 1.3·10 3 Ns m −3 , quadratic damping coefficient b 2 = 4.1 · 10 3 Nsm −2 and linear damping coefficient b 1 = 4.7 · 10 3 Ns m −1 . The matrix form of the homogeneous governing differential equation is

Applied Mechanics and Materials Vol. 885
[ In this mathematical model, the assumptions about the stiffness behavior in model b) and the damping behavior in model c) are combined. The resulting governing relation in the matrix form is

Solution of the Relative Compression Response of the 2DOF Mathematical Models a) to d)
The solution in time domain for (7), (10), (13) and (14) is described in a form that leads to the dynamic response y m,max in (3). It is to be noted that y m,max is applicable for the different linear, nonlinear, axiomatic and empiric assumptions for stiffness and damping behaviour. The solutions are obtained by numerical integration via R K algorithm. From (5) and (6), the accelerationsz u andz l arez which leads to the state-space form of (15) and (16)

Uncertainty in Mechanical Engineering III
Similarly, the general time domain solution to the governing relations for linear, nonlinear, axiomatic and empiric assumptions in (7), (10), (13) and (14) can be expressed in a generalized nonlinear state space form byẋ The summation of internal forces acting on m u and m l in (18) are described as internal force functions f u (x 1 , x 2 , x 3 , x 4 ) and f l ( x 1 , x 2 , x 3 , x 4 ), where x 1 , ..., x 4 are the state variables that represent the displacements and the velocities in (17). The functions f u (x 1 , x 2 , x 3 , x 4 ) and f l (x 1 , x 2 , x 3 , x 4 ) for the mathematical models a) to d) are listed in Tab. 1. The state vectorẋ is numerically integrated using a fourth order R K algorithm to obtain the desired state variables x 1 , x 2 , x 3 and x 4 . The numerical computation is started with the initial conditions in (1) and (2) along with the numerical time integration step t s = 0.001 s. Hence, the general desired dynamic output becomes Tab. 2 summarizes the model parameters for the mathematical models a) to d). According to Fig. 1 Applied Mechanics and Materials Vol. 885 9 with multi variable inputs such as the experimental mass of the upper structure m u and the experimental drop height h. The measured dynamic output variable in (20) comprises of the maximum compression z r,e,max , maximum spring damper force F sd,e,max and the maximum ground impact force F ef,e,max . Fig. 6 illustrates the MAFDS and the experimental test rig perspectively.   The experimental compression z r,e,max of the MAFDS is recorded by the displacement sensor 11 in Fig. 5 (a) that is based on the linear displacement variable transformer (LVDT) principle. The total load exerted F ef,e,max by the MAFDS on the ground is recorded by the three-axial strain gauge force sensor 12. Similarly, the force F sd,e,max exerted by the spring damper component of the MAFDS is recorded by the 1-axial strain gauge force sensor 13.

Uncertainty in Mechanical Engineering III
The experimental drop tests are designed using a simple full factorial (FF) design by varying the mass of the upper structure m u and the drop height h of the experimental MAFDS, [29]. The FF design takes into account the combinations within a certain range and amount of tests of independent variables and the order of experimental test of these combinations are randomized to avoid systematic errors. In this approach, the significance of the two independent variables m u and h, each having three different values, is considered, see Tab. 4. Fig. 7 shows the 3x3 FF design of experiments using m u and h as varying input variables. This design takes into account both the main effects and the interaction effects of the two variables m u and h. The main effect considers the influence of one independent variable on the output vector y e,max , while ignoring the influence of other independent variables. In Fig. 7, the experimental mass m u varies as per the combinations of (m u,1 ; h 1 ), (m u,2 ; h 1 ) and (m u,3 ; h 1 ), while the experimental drop height h remains constant. In this case, the variation in the experimental output y e,max is influenced by single variable m u . Conversely, an interaction effect of the independent variable occurs when its influence on the output relies on the value of the other independent variable also. For example, in the combinations (m u,2 ; h 3 ) and (m u,3 ; h 2 ) the output y e,max is influenced by both m u and h.

Quantification of Model Uncertainty
The four different simplified 2DOF mathematical models a) to d) introduced in section predict the dynamic response differently. The resulting model uncertainty is quantified with the B intervalbased hypothesis testing approach [20], [30]. This approach includes the uncertainty in the mathematical model prediction due to uncertain model parameters and the uncertainty in the experimental observations due to measurement variability.

Bayesian Interval Hypothesis-Based Method
This chapter gives an overview of the Bayesian interval-based approach to quantify the model uncertainty in the mathematical models. The first step is to estimate the difference data matrix Y d between the model simulation data matrix Y m,max and the experimental data matrix Y e,max . The difference data matrix Y d consists of difference vectors y d ∈ R 1×P , which is tested under the Bayesian interval based hypothesis approach. The number of columns P of the difference data matrix Y d represents the number of output variables which is equal to 3, according to (3). The null hypothesis H 0 occurs when the difference vector y d falls within the acceptable threshold limits ε. Conversely, the alternate hypothesis occurs when the difference vector y d falls outside the acceptable threshold limits ε. The Bayes factor B is computed by estimating the probability of the occurrence of the null hypothesis p(y d |H 0 ) to the probability of occurrence of the alternate hypothesis p(y d |H 1 ). Fig. 8 summarises the Bayesian interval approach to, first, derive the random difference row vector y d .  (21) The next step is to estimate the difference matrix Y d ∈ R D×P between the simulation data matrix Y m,max and the experimental data matrix Y e,max . The rows of the difference matrix Y d consist of the difference vector y d ∈ R 1×P . It is computed by taking the row vector difference the between y m,max and y e,max . Hence, the number of rows of Y d is obtained as D = M · E. The difference vector y d ∈ R 1×P in (22) is tested under the B interval based hypothesis approach for d = 1, 2, 3, ..., D.

Random Difference Row Vector
where the difference index The difference vector y d is assumed to follow a discrete Gaussian distribution with a mean vector µ d ∈ R 1×P and a covariance matrix Σ d ∈ R P ×P . The reason for assuming a Gaussian distribution is due to the fact that also the randomly picked samples y m,max and y e,max for the simulation and the experiment are assumed to follow a Gaussian distribution in these first investigations. However, sampling according to other distributions like e. g. uniform, Weibull or exponential distributions could be an option and are subject to further investigations.

Difference Row Vector under Null or Alternate Hypothesis
In section , the non-deterministic difference vector y d is defined in (22). Under the framework of the Bayesian interval hypothesis approach, the null hypothesis H 0 is represented as H 0 : |y d | ≤ ε, where ε is the threshold vector [20], [30]. The threshold vector ε is defined as where µ e is the mean of the experimental output vector y e and b th is the threshold coefficient that is to be defined by an expert or a model user [9]. Conversely, the alternate hypothesis is defined as H 1 : |y d | > ε. The subsequent subsection defines the likelihoods of observing the null hypothesis H 0 and the alternate hypothesis H 1 . The likelihood of null hypothesis is the probability of the non-deterministic difference vector y d falling within the the threshold limits [−ε, ε]. The covariance matrix Σ d ∈ R P ×P is computed from the difference matrix Y d . Conversely, the likelihood of the alternate hypothesis H 1 is

14
Uncertainty in Mechanical Engineering III

Estimation of Bayes factor B
According to the Bayes theorem, the bayes factor is the ratio of the the difference vector y d falling under the null hypothesis H 0 and the alternate hypothesis H 1 [31]. The Bayes factor is dependent on the prior assumptions of the model parameters k, b and k ef as well as on the threshold interval [−ε, ε].
As per the authors in [12], the model confidence M based on the Bayes factor B is determined as

Model Confidence for Models a) to d)
In this section the effect of the threshold vector ε (24) on the model confidence M (28) is discussed. The Bayes factor B (27) as well as its corresponding model confidence M are computed for all the four mathematical models a) to d). Fig. 9 illustrates the effect of variation of the threshold coefficient b th on the variation of the model confidence M for all the mathematical models a) to d). It is seen that the model b) with nonlinear stiffness and linear viscous damping shows lower model confidence compared to the models a), c) and d). For example, model b) has M = 80% for a threshold of around 47 %, whereas models a), c) and d) approach 80 % confidence for lower threshold values. The models c) and d) show similar variation of M with respect to the increasing value of b th . This implies that the nonlinear damping assumption in both cases c) and d) seems more adequate than the linear assumptions in a) and b) -regardless of linear or nonlinear stiffness assumption, and leads to higher confidence. Using (28), the model confidence approaches a value of 100 % when the Bayes factor B approaches ∞. Conversely, for a value of B = 1, the corresponding M = 50%.

Conclusion and Outlook
In this contribution, model uncertainty in a two degree of freedom mathematical model that is represented using four different functional relations is quantified via B interval based hypothesis approach. The four different models are based on linear, nonlinear, axiomatic and empiric assumptions. The B factor and its corresponding model confidence are used to quantify the model uncertainty. It is found that the mathematical models with empiric nonlinear damping show higher model confidence compared to the models with linear damping assumptions, regardless of linear or nonlinear stiffness assumption. In future work, the authors explore the possibility of including the uncertainty about the functional relations of the simple two degree of freedom mathematical model vs. a more complex finite element model under the framework of M C M C simulation.