Linear and Nonlinear Models for Drop Simulation of an Aircraft Landing Gear System with MR Dampers

Abstract

In this study, our focus is on the drop test simulation of an MR (Magnetorheological) damper-based main landing gear (MRMLG), aiming to explore multi-degree-of-freedom (DOF) dynamic models during aircraft landing. Three different 6-DOF dynamic models are proposed in this work, and their drop performances are compared with results achieved by commercial software. The proposed models include a nonlinear aircraft model (NLAM); a linearized approximated aircraft model (LAAM) linearizing from the nonlinear equations of motion in NLAM; and a fully approximated aircraft model (FAAM) which linearizes the MRMLG’s strut force model. In order to evaluate the drop performance of the aircraft landing gear system with MR dampers, a 7-DOF aircraft model incorporating the nonlinear MRMLG was formulated using RecurDyn. The principal comparative parameters are the coefficient of determination (R²) for the system response of each model with the RecurDyn model and root mean square error (RMSE), which is the ensemble of CG displacement data for each model. In addition, the ensemble of time series data is created for diverse drop scenarios, providing valuable insights into the performance of the proposed drop test models of an aircraft landing gear system featuring MR dampers.

1. Introduction

Drop simulation and testing play a key role in the design and evaluation of aerospace landing systems [1]. Accurately predicting the behavior of landing gear under various conditions is essential for ensuring the safety and reliability of aircraft operations [2]. Numerical modeling using computational tools has become an indispensable process in this domain, offering a cost-effective and efficient means of studying complex dynamic systems [3,4]. As aircraft technology continues to advance, there is a growing demand for more sophisticated landing gear systems that can adapt to varying landing conditions and provide precise attitude control. One promising approach is the use of magnetorheological (MR) dampers in the main landing gear to enhance its performance [5,6].

The MR damper-based main landing gear (MRMLG) system represents a significant leap forward in landing gear technology. By incorporating MR dampers, this innovative system enables the real-time adjustment of damping characteristics, providing precise control over landing dynamics. This adaptive characteristic allows for enhanced shock absorption, reduced vibrations, and improved stability during the critical phases of landing with sophisticated control strategies [7,8,9,10]. In the past 15 years, various attempts have been made to replace traditional landing gear with the MRMLG. For instance, a compact MRMLG was manufactured and pressure losses were compared between the derived model and experimental data under low-speed conditions [11,12]. Also, several technologies utilizing MRMLG have been developed, employing dimensionless parameters like the Bingham number to achieve efficient shock energy dissipation and ensure a smooth landing experience with consistent stroke regardless of the applied impulse [13,14]. In another study, the Bingham-plastic model was employed to model MRMLG for a helicopter, and an impact test was conducted using synthesized MR fluid, achieving a peak stroking load of up to 10 kN at a velocity of 4 m/s [15]. By eliminating the pneumatic chamber, a helicopter MRMLG design achieved a damping force of almost up to 40 kN, and the impact was effectively reduced through the application of various controls [16]. However, there are currently few studies that have derived an entire aircraft model by incorporating MRMLG and conducted a comparison and a validation of the model using a commercial dynamics analysis program. Moreover, to fully exploit the potential of MRMLG systems, accurate dynamic models are required to capture the complex interactions between the aircraft, landing gear, and control systems. In other words, modeling a multi-degree-of-freedom (DOF) aircraft equipped with the MRMLG is significant, as it provides a comprehensive understanding of the landing behavior and performance. By accurately representing the behavior of the MRMLG system, the model can provide valuable insights into its dynamics, shock absorption capabilities, and overall landing performance. This modeling approach enables engineers and researchers to evaluate the effectiveness of the MRMLG system, optimize its design parameters, and explore potential enhancements for improved aircraft stability and passenger comfort during landing operations. In the field of multi-DOF aircraft modeling with the MRMLG, there have been some notable works. One study proposed the landing gear’s 6-DOF vibrational model and state-space equation, followed by the design of a control system comprising admittance and displacement control loops [17]. In another study, a 3-DOF aircraft model was developed to capture the landing impact’s bounce, pitch, and roll motions. The controlled responses were analyzed comparatively using the robust H_∞ and LQR control strategies [18]. However, it should be noted that both papers suffer from the limitations of model accuracy and reliability for practical application, as they employed simplistic linear models for the aircraft and MR damper.

This work aimed to develop and validate various 6-DOF aircraft models to evaluate precise dynamic drop motions during aircraft landings. We propose three different models: the nonlinear aircraft model (NLAM), the linear approximated aircraft model (LAAM), and the fully approximated aircraft model (FAAM). These models are constructed based on the nonlinear equations of motion and linear approximations, considering the effect of MR dampers on the strut force. To validate the accuracy of the proposed models, RecurDyn, which is known as a powerful multi-body analysis tool, was utilized. A comprehensive comparison was made between the models and the RecurDyn model, considering evaluation metrics such as the coefficient of determination (R²) and the root mean square error (RMSE). Additionally, an ensemble of time series data was created to analyze the performances of the models under various drop scenarios.

The outcomes of this research will contribute to the development of effective control strategies for MRMLG, and will also enhance our understanding of the behavior of landing gear during drop motions. The validation and analysis of the NLAM, LAAM, and FAAM models will assist in selecting an appropriate model for dynamic simulations which can provide accurate predictions of MRMLG performance. The findings in this work will be valuable for improving the design and implementation of MRMLG systems to achieve enhanced aircraft landing performances. This research will also serve as a foundation for further advancements in modeling techniques and their applications in aerospace engineering.

2. MRMLG Design and Drop Test

The MRMLG consists of the structure shown in Figure 1. The chambers consist of three sections: Chambers 1 and 2, as well as the orifice, are filled with MR fluid, and the gas chamber is mainly filled with nitrogen gas. MR fluid can flow freely into chambers 1 and 2 through the orifice. The biggest feature of the MRMLG is that the flow resistance of the fluid can be controlled by the electromagnet around the orifice. The electromagnet consists of inner and outer cores. The solenoid is wound around the groove of the inner core, and the orifice beside the electromagnet should be located in the fully developed region [6]. Therefore, the entry is installed on the top of the electromagnet to ensure that the orifice corresponds to the entrance region. The electromagnetic is composed of multiple cores, and it allows for a more efficient magnetic flux density. Chamber 2 and the gas chamber are separated by the separator. The gas in the gas chamber serves as a pneumatic spring that converts kinetic energy into elastic energy when the aircraft is landing, creating a repulsive force depending on the stroke of the rod. Seals divide the chamber by suppressing the penetration of fluid, and bushes mainly guarantee the linear motion of each part. The trunnion connects to the fuselage’s gearbox. During landing, most of the impact force acting on the runway transfers to the wheels and rods, and the MR fluid dampens the force acting on the main strut. The braces and columns dissipate the forces acting on the fuselage so that the aircraft can land safely.

The strut force model of the MRMLG was derived in detail from previous studies, and the derived equation is as follows [6]:

$$F_s=F_{gas}(s)+F_{hyd}(\dot{s})+F_{mr}(\dot{s},\tau )+F_c$$

(1)

where Fs is the strut force, and F_gas, F_hyd, F_MR, and F_c are the pneumatic force, hydraulic force, field-dependent force, and Coulomb friction force, respectively. In this study, the Coulomb friction was ignored, since it is very small compared to other forces. Each force can be calculated as:

$$F_{gas}(s)=\left(P_{atm}+P_0\right)A_2{\left(\frac{V_0}{V_0-A_{2}s}\right)}^{\gamma }-P_{atm}{A}_2$$

(2)

$$F_{hyd}(\dot{s})=\frac{12\eta l{A}_1}{b{d}^3}\dot{s}+K\frac{\rho A_1{}^2}{2b^2{d}^2}{\dot{s}}^2$$

(3)

$$F_{mr}(\dot{s},\tau )=\left(2.07+\frac{30\eta A_1\left|\dot{s}\right|}{30\eta A_1\left|\dot{s}\right|+b{d}^2\left|\tau \right|}\right)\frac{l_p}{d}\left|\tau \right|\tanh \left(\frac{\dot{s}}{\epsilon }\right)$$

(4)

where P_atm is the atmospheric pressure; P₀ and V₀ are the initial gas pressure and volume; γ is the ratio of the specific heats; A₁ and A₂ are the inner diameter of the main strut and the outer diameter of the rod, respectively; η and ρ are the dynamic viscosity and density of the MR fluid; L, b, and d are the length, circumference, and gap, respectively; and L_p is the total length of the electromagnet, except the solenoid length. K is the total sum of system losses caused by the turbulence of fluid around the orifice. K was calculated depending on the flow path gap size using computational fluid dynamics in a previous work [19]. ε is a small, positive value that is used to adjust the rate of the hyperbolic tangent function, tanh. $s$ and $\dot{s}$ are the stroke and the velocity of the rod, respectively. τ is the yield stress, which is a function of the magnetic field intensity, and the yield stress characteristics are provided by the MR fluid manufacturer [20].

The MRMLG can be modeled as a two-degree-of-freedom (2-DOF) by imitating the drop test conditions. In general, half of the weight of the fuselage is set as sprung mass, m_s, and the weight of the rod, electromagnet, tire, and wheel can be replaced by unsprung mass, m_u. The equations of motion for sprung and unsprung masses are:

$${\ddot{z}}_s=g-F_s/m_s$$

(5)

$${\ddot{z}}_u=g+F_s/m_u-F_t/m_u$$

(6)

where g is the gravitational acceleration and ${\ddot{z}}_s$ and ${\ddot{z}}_u$ are the sprung and unsprung accelerations, respectively. F_t can be calculated in the tire force as:

$$F_t=\left(k_t{z}_u-k_{t}r\right)\mathscr{H}(z_u-r)$$

(7)

where k_t and r are the tire stiffness and ground displacement, respectively; H is the Heaviside step function. The drop test is conducted under the condition of no ground displacement. To verify the model via the drop test, a drop simulation was performed numerically by building a 2-DOF model based on Equations (5) and (6). The parameters used in the drop simulation are listed in

Table 1. MRMLG parameters and physical quantities.

Parameter	Notation	Value	Unit
Atmospheric Pressure	Patm	1.013	bar
Initial Gauge Pressure	P0	4.000	bar
Initial Volume of the Gas Chamber	V0	454.0	cm3
Polytropic Index of the Nitrogen	γ	1.300	-
Inner Area of the Main Strut	A1	25.52	cm3
Outer Area of the Rod	A2	20.19	cm3
Viscosity of MRF-132DG [20]	η	112.0	mPa s
Density of MRF-132DG [20]	ρ	3.050	g/cm3
Orifice Length	l	130.0	mm
Orifice Perimeter	b	139.4	mm
Orifice Gap	d	1.300	mm
Effective Pole Length	lp	49.40	mm
Total Loss Coefficient	K	2.836	-
Slope Value of the Hyperbolic Tangent	ε	0.050	-
Sprung Mass	ms	680.0	kg
Unsprung Mass	mu	18.00	kg
Tire Stiffness	kt	412.0	kN/m
Gravitational Acceleration	g	9.807	m/s2

.

The function of magnetic field intensity for the input current was obtained in the previous study through the ANSYS MAXWELL electromagnetic analysis program [19]. Based on the determined parameters, the yield stress for the input current was calculated:

$$\tau =40.5\tanh {(1.3I)}^{1.8}$$

(8)

The range of I in Equation (8) is from 0 to 2 A and the units of τ and I are kPa and A, respectively.

A drop test was performed to validate the MRMLG 2-DOF model, as shown in Figure 2. The MRMLG and dummy mass hung on the drop carriage. The dummy mass was used to match the sprung mass to 680 kg. The MRMLG was equipped with a wire sensor and a laser sensor to monitor the stroke and sprung mass displacement, respectively, and pressure sensors were installed on chambers 1 and 2 to indirectly measure the strut force. The DAQ acquired sensor data and sent a signal to the power amplifier. The power amplifier then applied the input current to the MRMLG based on the signal. As a drop test condition, a current of 0 A to 2 A was applied to the MRMLG at a sink rate of 3.050 m/s just before the tire contacted the ground.

Figure 3 shows the simulation and test results of the strut force in relation to the stroke velocity. The solid line represents the simulation results, while the dots represent the test results. The measured strut force was obtained by using the pressure sensors attached to chambers 1 and 2. The accuracy of the measurements and simulations was evaluated using the root mean square error (RMSE) and coefficient of determination (R²). RMSE is the square root of the mean square error of two datasets, and R² measures the proportion of the variance in the dependent variable that can be explained by the independent variable, determined in the range of 0 to 1. As the value increases, the correlation between the two datasets increases. The hydraulic force model was relatively accurate when off-state, since the RMSE and R² were 1.018 kN and 0.950. the RMSE and R² of 1 A were 2.179 kN and 0.844, showing some model error at high speeds. The test results of 2 A were higher than the predicted values in all velocity ranges, which were measured with an RMSE of 2.309 and an R² of 0.746, respectively. The simulation and test results were similar, indicating that the modeling of MRMLG had been conducted correctly.

Figure 4 illustrates the spring-damping characteristics of MRMLG. Figure 4a shows the strut force exerted by the MRMLG as a function of stroke when a 0 A input current was applied. In the simulation, the maximum stroke and strut force were predicted to be 199.9 mm and 28.30 kN, respectively, while the actual measurements yielded 205.6 mm and 28.24 kN. The point of equilibrium in the simulation and test were determined to be 165.7 mm and 192.7 mm, respectively. The discrepancy in the equilibrium point can be attributed to the compression process of the gas chamber, which underwent an adiabatic process during the initial impact but gradually released heat to the surrounding areas, transitioning to an isothermal process over time. The RMSE values for the stroke and strut force were 29.31 mm and 1.556 kN, respectively, with corresponding R² values of 0.738 and 0.889. Figure 4b describes the characteristics of 2 A input. The predicted and measured values for strut force were 29.01 kN and 29.97 kN, respectively, with corresponding stroke values of 174.9 mm and 176.53 mm. The RMSE and R² values for stroke were 10.56 mm and 0.792, respectively, while for strut force they were 2.003 kN and 0.768, respectively. Irrespective of the input current, the strut force during the initial impact (within a stroke range of 0 to 200 mm) showed similarity between the model and test results. However, it has been verified that there was some model error for the strut force after reaching maximum stroke.

Figure 5 displays the predicted (solid line) and measured (dashed line) system responses to an impact. (a) and (b) show the sprung displacement (z_u) and unsprung displacement (z_s) for 0 A and 2 A current inputs, respectively. For the off-state, the average values of z_s and z_u were predicted to be 176.1 mm and 17.50 mm, while the measured values were 193.3 mm and 13.25 mm, respectively. The predicted final values were 184.9 mm and 17.00 mm, whereas the measured final values were 11.65 mm and 204.0 mm. The RMSE values were 21.10 mm and 6.667 mm for z_s and z_u, respectively, with corresponding R² values of 0.890 and 0.851. For the on-state, the average values of z_s were 183.1 mm and 170.6 mm in the simulation and the test, with equilibrium points of 184.5 mm and 182.6 mm, respectively. The RMSE and R² values for z_s were 18.55 mm and 0.701. The predicted and measured average values for z_u were 18.29 mm and 9.182 mm, with final values of 17.29 mm and 8.621 mm, respectively. The RMSE and R² values for z_u were 10.78 mm and 0.693. In the simulation, the absence of friction led to chattering in z_s and z_u, which was not observed in the actual measurements due to the presence of friction. As a result, the R² value was slightly lower than that of the off-state. The damping elements of the tire and structures were not included in the model, as observed in Equations (5) and (6). Consequently, there was some discrepancy between the model and test results after the initial stroke.

3. Nonlinear Model and Linear Approximation Model

For the analysis of the aircraft’s landing motion, previous research had derived the kinematics and dynamics of the 6-DOF aircraft landing system using the Euler–Lagrange equation [21]. Figure 6 shows the ground’s absolute coordinate system (x, y, z) and the lumped aircraft model’s relative coordinate system (u, v, w); position vectors; and roll (φ), pitch (θ), and yaw (ψ) directions [22]. The aircraft model’s generalized coordinate q and i-th position vector P_i were defined in the previous research:

$$q={\left[\begin{array}{cccccc}{z}_0& s_1& s_2& s_3& \phi & \theta \end{array}\right]}^{\mathrm{T}}$$

(9)

$$P_i={\left[\begin{array}{ccc}{X}_i& Y_i& Z_i\end{array}\right]}^{\mathrm{T}}$$

(10)

where s₁ represents the stroke of the nose gear; s₂ and s₃ are the strokes of the left and right main gears; and φ and θ are the roll and pitch angles, respectively. X_i, Y_i, and Z_i denote the distances along the i-th x-, y-, and z-axes, respectively. Ψ, shown in Figure 6, is the yaw angle, which arises from the steering of the nose gear. This study assumes no yaw motion and therefore ignores it. Moreover, the equation of motion for a nonlinear aircraft model (NLAM) is as follows:

$$\ddot{q}=M_G{}^{-1}\left(Q-\left(f_{ℒ}-M_G\ddot{q}\right)\right)$$

(11)

where M_G is the generalized mass matrix, Q is the generalized force vector, and f_L is the solution to the Euler–Lagrange equation. The equation of motion in Equation (11) has a disadvantage in that the formula is complicated and it is difficult to apply the linear control logic. If the formula is simplified by linear approximation, the system analysis and linear control logic become simpler, but the accuracy of the model decreases.

The nonlinear position vectors and the equation of motion, derived from kinematics and dynamics in the previous research, are linearly approximated. During linearization, it is assumed that the perturbations in the system variables are very small. Additionally, the sinusoidal function is also linearly approximated, since the perturbation of the angle is also very small. First, the position vector in Figure 6 can be approximated through the Jacobian matrix as follows:

The position vectors calculated through the linearization were:

$$p_0={\left[\begin{array}{ccc}0& 0& z_0\end{array}\right]}^{\mathrm{T}}$$

(12)

$$p_1=\left[\begin{array}{c}{U}_n+W_n\theta \\ -W_n\phi \\ W_n+z_0-U_n\theta \end{array}\right]$$

(13)

$$p_2=\left[\begin{array}{c}-U_m+W_m\theta \\ V_m-W_m\phi \\ W_m+z_0+V_m\phi +U_m\theta \end{array}\right]$$

(14)

$$p_3=\left[\begin{array}{c}-U_m+W_m\theta \\ -V_m-U_m\phi \\ W_m+z_0-V_m\phi +U_m\theta \end{array}\right]$$

(15)

$$p_4=\left[\begin{array}{c}{U}_n+L_n{\theta }_n-{\theta }_n{s}_1+\left(L_n+W_n\right)\theta \\ -\left(L_n+W_n\right)\phi \\ L_n+W_n+z_0-s_1-\left(U_n+L_n{\theta }_n\right)\theta \end{array}\right]$$

(16)

$$p_5=\left[\begin{array}{c}-U_m+L_m{\theta }_m-{\theta }_m{s}_2+\left(L_m+W_m\right)\theta \\ V_m-\left(L_m+W_m\right)\phi \\ L_m+W_m+z_0-s_2+V_m\phi +\left(U_m-L_m{\theta }_m\right)\theta \end{array}\right]$$

(17)

$$p_6=\left[\begin{array}{c}-U_m+L_m{\theta }_m-{\theta }_m{s}_3+\left(L_m+W_m\right)\theta \\ -V_m-\left(L_m+W_m\right)\phi \\ L_m+W_m+z_0-s_3-V_m\phi +\left(U_m-L_m{\theta }_m\right)\theta \end{array}\right]$$

(18)

where z₀ is the distance from the origin of the absolute coordinate system to the fuselage’s center of gravity (CG). Other geometric parameters used in the position vectors are listed in

Table 2. Geometric parameters for the benchmarked aircraft model.

Parameter	Notation	Value	Unit
y-angle of the nose gear	θn	12.8	degree
y-angle of the main gear	θm	6.5	degree
x-distance of the nose gear from the CG	Un	1.652	m
z-distance of the nose gear from the CG	Wn	0.124	m
x-distance of the main gear from the CG	Um	0.393	m
y-distance of the main gear from the CG	Vm	1.479	m
z-distance of the main gear from the CG	Wm	0.038	m
Total length of fully extended nose gear	Ln	0.756	m
Total length of fully extended main gear	Lm	0.830	m

.

Second, assuming that the nominal initial condition and perturbation of the generalized coordinates are very small and time-invariant, M_G can be expressed as:

$$\begin{array}{l}{M}_L={M_G|}_{q=0}=\\ \left[\begin{array}{cccccc}{m}_0+m_4+m_5+m_6& -m_4& -m_5& -m_6& \left(m_5-m_6\right)V_m& \begin{array}{l}-m_4\left(U_n+L_n{\theta }_n\right)\\ +\left(m_5+m_6\right)U_m\\ -\left(m_5+m_6\right)L_m{\theta }_m\end{array}\\ -m_4& m_4& 0& 0& 0& m_4\left(U_n+W_n{\theta }_n\right)\\ -m_5& 0& m_5& 0& -m_5{V}_m& -m_5\left(U_m+W_m{\theta }_m\right)\\ -m_6& 0& 0& m_6& m_6{V}_m& -m_6\left(U_m+W_m{\theta }_m\right)\\ \left(m_5-m_6\right)V_m& 0& -m_5{V}_m& m_6{V}_m& \begin{array}{l} I_x+m_4{\left(L_n+W_n\right)}^2\\ +\left(m_5+m_6\right)V_m^2\\ +\left(m_5+m_6\right){\left(L_m+W_m\right)}^2\end{array}& \begin{array}{l} \left(m_5-m_6\right)U_m{V}_m\\ -\left(m_5-m_6\right)L_m{V}_m{\theta }_m\end{array}\\ \begin{array}{l}-m_4\left(U_n+L_n{\theta }_n\right)\\ +\left(m_5+m_6\right)U_m\\ -\left(m_5+m_6\right)L_m{\theta }_m\end{array}& \begin{array}{l} m_4{U}_n\\ -m_4{W}_n{\theta }_n\end{array}& \begin{array}{l}-m_5{U}_m\\ -m_5{W}_m{\theta }_m\end{array}& \begin{array}{l}-m_6{U}_m\\ -m_6{W}_m{\theta }_m\end{array}& \begin{array}{l} \left(m_5-m_6\right)U_m{V}_m\\ -\left(m_5-m_6\right)L_m{V}_m{\theta }_m\end{array}& \begin{array}{l} I_y+m_4{U}_n^2\\ +m_4{\left(L_n+W_n\right)}^2\\ +2m_4{L}_n{U}_n{\theta }_n\\ +\left(m_5+m_6\right)U_m^2\\ +\left(m_5+m_6\right){\left(L_m+W_m\right)}^2\\ -2\left(m_5+m_6\right)L_m{U}_m{\theta }_m\end{array}\end{array}\right]\end{array}$$

(19)

Third, the vectors of the equation of motion in Equation (11) can be linearly approximated through the Jacobian matrix as follows [23]:

$$\begin{array}{l}{Q}_L={J_Q(q)|}_{\begin{array}{l}q=0\\ F=0\end{array}}q+{J_Q(F)|}_{\begin{array}{l}q=0\\ F=0\end{array}}F\\ =\left[\begin{array}{cccccc}0& 0& 0& -1& -1& -1\\ -1& 0& 0& 1& 0& 0\\ 0& -1& 0& 0& 1& 0\\ 0& 0& -1& 0& 0& 1\\ 0& 0& 0& 0& -V_m& -V_m\\ 0& 0& 0& U_n+L_n{\theta }_n& -U_m+L_m{\theta }_m& -U_m+L_m{\theta }_m\end{array}\right]F\end{array}$$

(20)

$$\begin{array}{l}{g}_L/g=\left({J_{\left(f_{ℒ}-M_G\ddot{q}\right)}(q)|}_{\begin{array}{l}q=0\\ \dot{q}=0\end{array}}q+{J_{\left(f_{ℒ}-M_G\ddot{q}\right)}(\dot{q})|}_{\begin{array}{l}q=0\\ \dot{q}=0\end{array}}\dot{q}+{\left(f_{ℒ}-M_G\ddot{q}\right)|}_{\begin{array}{l}q=0\\ \dot{q}=0\end{array}}\right)/g\\ = \left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -m_4{\theta }_n\\ 0& 0& 0& 0& 0& -m_5{\theta }_m\\ 0& 0& 0& 0& 0& -m_6{\theta }_m\\ 0& 0& 0& 0& \begin{array}{l} m_4\left(L_n+W_n\right)\\ +\left(m_5+m_6\right)\left(L_m+W_m\right)\end{array}& 0\\ 0& -m_4{\theta }_n& -m_5{\theta }_m& -m_6{\theta }_m& 0& \begin{array}{l} m_4\left(L_n+W_n\right)\\ +\left(m_5+m_6\right)\left(L_m+W_m\right)\end{array}\end{array}\right]q\\ +\left[\begin{array}{c}-m_0-m_4-m_5-m_6\\ m_4\\ m_5\\ m_6\\ -\left(m_5-m_6\right)V_m\\ m_4\left(U_n+L_n{\theta }_n\right)-\left(m_5+m_6\right)\left(U_m-L_m{\theta }_m\right)\end{array}\right]\end{array}$$

(21)

where F is the force vector, i.e.,

$$F={\left[\begin{array}{cccccc}{F}_{s,1}& F_{s,2}& F_{s,3}& F_{t,4}& F_{t,5}& F_{t,6}\end{array}\right]}^{\mathrm{T}}$$

(22)

The elements of the force vector were arranged as follows according to Equations (1)–(4) and (7)–(8), respectively:

$$F_{s,1}=k_1{s}_1+c_1{\dot{s}}_1$$

(23)

$$F_{s,2}=F_{gas}(s_2)+F_{hyd}({\dot{s}}_2)+F_{mr}({\dot{s}}_2,I_2)$$

(24)

$$F_{s,3}=F_{gas}(s_3)+F_{hyd}({\dot{s}}_3)+F_{mr}({\dot{s}}_3,I_3)$$

(25)

$$F_{t,4}=k_4\left(z_4-r_4\right)\mathscr{H}(z_4-r_4)$$

(26)

$$F_{t,5}=k_5\left(z_5-r_5\right)\mathscr{H}(z_5-r_5)$$

(27)

$$F_{t,6}=k_6\left(z_6-r_6\right)\mathscr{H}(z_6-r_6)$$

(28)

where k₁ and c₁ represent the spring coefficient and damping coefficient of the nose gear strut; I₂ and I₃ denote the control inputs (current) applied to the right and left main gears; k₄, k₅, and k₆ represent the tire stiffness of the nose, right, and left gears; and r₄, r₅, and r₆ refer to the runway roughness transmitted to the tires of the nose, right, and left gears, respectively. In this study, it was assumed that the aircraft would be landing on an ideal runway, so the runway’s roughness was neglected. As a result, the motion equation of a linear approximated aircraft model (LAAM) was:

$$\ddot{q}=M_L^{-1}\left(Q_L-g_L\right)$$

(29)

Although the affine equations of motion for the aircraft landing system were derived, the pneumatic spring, hydraulic force, and field-dependent force of the MRMLG in Equations (2)–(4) have not yet been linearized. As a result, the system equation of the LAAM satisfied the linear approximation, but the elements of the generalized force vector included the nonlinear MRMLG model. Therefore, by linearizing the MRMLG, it is possible to construct a fully approximated aircraft model (FAAM). For the strut force model in Equations (2)–(4), a linearized relationship can be obtained using the linear regression technique, as shown in Figure 7.

In Figure 7, k, c, and h are the linearized pneumatic spring coefficient, damping coefficient, and coefficient for the input current of the MRMLG, respectively. The following linear regressed forces were used instead of Equations (24) and (25) in the FAAM model:

$$F_{s,2}=k_2{s}_2+c_2{\dot{s}}_2+h_2{I}_2$$

(30)

$$F_{s,3}=k_3{s}_3+c_3{\dot{s}}_3+h_3{I}_3$$

(31)

All physical quantities used in NLAM, LAAM, and FAAM, as well as their descriptions, are summarized in

Table 3. Physical parameters for the benchmarked aircraft model.

Parameter	Notation	Value	Unit
Sprung (fuselage) mass	m0	1633	kg
Unsprung (nose and main gears) mass	m4, m5, m6	18.00	kg
Moments of inertia about the x-axis	Ix	6132	kg m2
Moments of inertia about the y-axis	Iy	2555	kg m2
Moments of inertia about the z-axis	Iz	8595	kg m2
Stiffness of the nose gear	k1	30.00	kN/m
Damping coefficient of the nose gear	c1	3000	N s/m
Linearized stiffness of the main gears	k2, k3	46.68	kN/m
Linearized damping coefficient of the main gears	c2, c3	7337	N s/m
Linearized coefficient of the control input	h2, h3	7439	N/A
Tire stiffness of the nose gear	k4	200.0	kN/m
Tire stiffness of the main gears	k5, k6	412.0	kN/m

.

4. Drop Simulation Comparison with RecurDyn Model

An MRMLG-based aircraft drop simulation was performed via the RecurDyn program to analyze the accuracy of the model. RecurDyn is a program that specializes in multi-body dynamics analysis. As shown in Figure 8, the aircraft model consisted of the fuselage, two MRMLGs, and one nose landing gear; each landing gear was composed of a rod, cylinder, and tire. The constraint settings were configured as follows: the fuselage, cylinder, rod, and tire were fixed; the cylinder and rod were allowed translation; and the tire and ground had a frictionless planar condition. The forces acting on the tire were set to act vertically upward from the ground. The friction force in the translational condition was neglected. Damping forces, pneumatic forces, and reaction forces due to tire were all applied based on the mathematical model in Equations (2)–(4). The required state variables of 7-DOF for modeling were obtained directly from the 3D model, including the CG of the aircraft and the landing gear strokes (nose, right, left gears), as well as the roll, pitch, and yaw of the landing gear.

Figure 9a,b illustrate the aircraft landing behavior during the free fall simulations performed using RecurDyn and MATLAB. The simulation settings were configured with a time step of 1 ms for both RecurDyn and MATLAB, and data were recorded for a total duration of 4 s. The simulation solver in RecurDyn was set to dynamic analysis, and the specific parameters which were used are summarized in Table 2 and Table 3. The NLAM, LAAM, and FAAM models were numerically analyzed using MATLAB’s ode45 solver. For the initial conditions, the sink rate was fixed at 3.050 m/s, while the roll and pitch angles and the input current were adjusted. The roll angles were set to 0°, 3°, and 6°, and the pitch angles were set to 0°, 5°, and 10°. Input currents of 0 A, 1 A, and 2 A were applied. In this study, the notation RφPθuI (pitch angle: φ degrees, roll angle: θ degrees, input current: I A) was used to define the initial conditions. For example, the notation R6P10u2 represents the initial condition with roll and pitch angles of 6° and 10°, respectively, and an input current of 2A. During the analysis, system responses and strut force characteristics were compared and analyzed for different models, including R0P0u0, R0P0u2, R6P10u0, and R6P10u2 of initial conditions. The reference model was RecurDyn, and the accuracy of the NLAM, LAAM, and FAAM models was evaluated based on R². Subsequently, ensemble processing was performed on the system responses for various initial conditions to evaluate the models. R² and RMSE were calculated with respect to RecurDyn to assess the models’ accuracy levels and errors.

Figure 10 shows the system’s response to the shock under the R0P0u0 condition. Figure 9a represents the displacement of the CG, and Figure 9b shows the pitch angular displacement. The roll angle remained constant at zero for all models. The R² between the RecurDyn model and the NLAM model was calculated as 1.000 for both the CG displacement and the pitch angle, indicating no difference. The R² values between the RecurDyn model and the LAAM model were 0.999 and 0.973, respectively, for the CG displacement and the pitch angle. However, the FAAM model showed very low R² values of 0.119 and 0.014, indicating a lack of correlation with the reference model. Figure 11 illustrates the system response under the R0P0u2 condition. In Figure 10a, the CG displacement, which yielded R² values of 0.998, 0.996, and 0.482 for the NLAM, LAAM, and FAAM models, respectively, is shown. In Figure 10b, the pitch angle is shown to exhibit R² values of 1.000, 0.973, and 0.941, in the same order. These results indicate an improvement in the accuracy of the FAAM model compared to the R0P0u0 condition.

Figure 12 depicts the strut force characteristics for the stroke under the R0P0u0 and R0P0u2 conditions. It is worth noting that the results exhibit similarities to the 2-DOF drop test simulation and test results shown in Figure 4. In Figure 12a, the maximum stroke values for RecurDyn, NLAM, LAAM, and FAAM are shown as 199.8 mm, 199.8 mm, 199.9 mm, and 252.3 mm, respectively, and the corresponding maximum strut force values as 27.66 kN, 27.69 kN, 27.97 kN, and 21.98 kN. In the figure, FAAM closely matches the other models up to a stroke of 180 mm, but diverges in accuracy beyond that point. This discrepancy was due to a model error in the pneumatic model for strokes above 180 mm, as depicted in Figure 7a. In the time series analysis, the R² values for stroke, in the order of NLAM, LAAM, and FAAM, were 1.000, 0.998, and 0.050, respectively. For strut force, the R² values were 1.000, 0.998, and 0.648, respectively. These results indicate that the accuracy of the FAAM model was significantly lower compared to NLAM and LAAM. In Figure 12b, under the R0P0u2 condition, the maximum stroke and strut force were approximately 176 mm and 29 kN for NLAM and LAAM, while for FAAM they were 140.3 mm and 32.35 kN, respectively. The deviation of approximately 3 kN in the maximum strut force for FAAM is attributed to an error in the input current model shown in Figure 7c, which consequently led to a shorter maximum stroke of approximately 36 mm compared to NLAM and LAAM. The R² values for stroke and strut force in the time series analysis were 1.000 and 0.996 for NLAM, 1.000 and 0.991 for LAAM, and 0.958 and 0.587 for FAAM, indicating an overall improvement in the accuracy of the FAAM model with the increasing input current.

The simulation results are compared between the 2-DOF model in Figure 7 and the R0P0u0 condition of the 7-DOF model using RecurDyn in Figure 12. Under 0 A, the R² value and RMSE for the stroke were 0.890 and 11.55 mm, respectively, while those for the strut force were 0.990 and 0.518 kN, respectively. Under 2 A, the R² values for the stroke and strut force were 0.999 and 0.727, respectively, while the RMSE for stroke was 1.334 mm and that for strut force was 2.440 kN.

Figure 13 and Figure 14 present the model responses and strut force–stroke characteristics for the R6P10u0, respectively, since the roll angle being non-zero leads to different responses between the right and left main gears in Figure 14. Compared to the R0P0u0, the NLAM model’s accuracy under the R6P10u0 condition was not significantly different, with R² values above 0.98 for all variables in the time series analysis. However, the LAAM model showed a notable decrease in accuracy. FAAM exhibited unreliable model accuracy.

Table 4. R² for system response of each model with the RecurDyn model (R6P10u0).

Variables	NLAM	LAAM	FAAM
Displacement of CG	0.998	0.876	0.086
Roll angle	0.999	0.846	0.299
Pitch angle	0.998	0.766	0.307
Stroke of the right gear	0.994	0.686	0.116
Stroke of the left gear	0.996	0.634	0.225
Strut force of the right gear	0.987	0.510	0.486
Strut force of the left gear	0.991	0.484	0.445

summarizes the R² calculations for the time series responses of each model under the R6P10u0 condition. In Figure 13a–c, a data-shifting phenomenon can be observed in LAAM, resulting in R² values below 0.9. The strut force of the right main gear had the lowest R² value of 0.510, which can also be attributed to the data shifting between the LAAM model and the time series data from RecurDyn. The R² values for the time series responses in FAAM do not exceed 0.5, with the stroke of the right main gear being at its lowest at 0.116.

Figure 15 and Figure 16 illustrate the CG displacements, angular displacements, and left–right characteristics of the main landing gear for the R6P10u2 condition. The R² values for the time series data are summarized in

Table 5. R² for system response of each model, along with the RecurDyn model (R6P10u2).

Variables	NLAM	LAAM	FAAM
Displacement of CG	0.988	0.960	0.681
Roll angle	0.992	0.939	0.509
Pitch angle	0.998	0.838	0.629
Stroke of the right gear	0.994	0.929	0.790
Stroke of the left gear	0.998	0.989	0.961
Strut force of the right gear	0.932	0.718	0.103
Strut force of the left gear	0.861	0.621	0.069

. In the NLAM model, the lowest R² value in the time series data for the R6P10u2 condition was 0.861 for the force of the left main strut, slightly lower compared to the R0P0u2 condition. NLAM is modeled with 6-DOF, while RecurDyn utilizes a 7-DOF model including yaw. It is observed that as the roll and pitch angles increased, the yaw angle fluctuated within 0.5 degrees. This fact suggests that the increase in model error may have contributed to the decrease in R² for the NLAM model. The R² values for the LAAM model range from 0.621 to 0.989 for the time series data, and the data-shifting phenomenon observed in R6P10u0 occurs only in Figure 15b. For the FAAM model, the R² values for the time series data were above 0.5 for all variables except the strut force, and data-shifting can be observed in Figure 16a–c. The minimum R² was calculated as 0.069 for the strut force. The reason why the strut force of the FAAM model overlaps with other models below 180 mm in Figure 11 is relatively straightforward. Figure 6a presents the pneumatic force model and linear regression results for stroke, showing a relatively consistent difference between the results up to approximately 180 mm. As a result, in Figure 11, Figure 13 and Figure 15, we can confirm that there is almost no error between FAAM and other models within 180 mm. In other words, the FAAM model closely overlaps with the other models.

To analyze the synthesized R² values of 4 s time series variables for various initial conditions, an ensemble of time series data was created for 27 initial conditions, resulting in a total of 108 s of data. Figure 17 displays the ensembles of CG displacement data for each model. The ensembles were used to perform evaluations of the coefficient of determination and to calculate the RMSE. Figure 18 shows the results of calculating the R² values by constructing ensembles for each time series variable and model, displayed in a bar chart format. The NLAM model showed R² values of 0.93 or higher for all variables, which is comparable to the RecurDyn model. The LAAM exhibited displacement responses with R² values above 0.85 and force models with R² values above 0.7, indicating its suitability for the drop simulation. However, the FAAM model showed R² values below 0.45 for all variables except hydraulic force. These results indicate that the FAAM is unsuitable for use as a drop simulation model.

Figure 19 displays histograms analyzing the model-specific RMSE for each variable in the ensemble. It can be observed that models with high R² values in Figure 18 exhibited relatively low RMSE values. Figure 19a displays the RMSE for strokes and CG displacement. The model-specific RMSE values were within 2 mm for NLAM and within 12 mm for LAAM, and reached a maximum of approximately 40 mm for FAAM, indicating the highest error range. Figure 19b shows the RMSE values for the roll and pitch angles. The RMSE for angular displacement was less than 0.1 degrees for NLAM and less than 1 degree for LAAM, but exceeded 1.5 degrees for FAAM. Figure 19c presents the RMSE values for the pneumatic and hydraulic forces of the right and left main gears, as well as field-dependent and strut forces. The RMSE for hydraulic force was below 1.5 kN for all models. The RMSE for pneumatic force was less than 0.4 kN for NLAM, less than 2.3 kN for LAAM, and below 4 kN for FAAM. The field-dependent and strut forces had RMSE values below 0.8 kN for NLAM and below 1.5 kN and 2.8 kN for LAAM, respectively. FAAM exhibited the highest RMSE values, approximately 3 kN for field-dependent force and reaching a maximum of 4.691 kN for strut force, indicating the largest degree of error.

5. Conclusions

This study focused on the dynamic drop simulation models of an aircraft landing system equipped with MR dampers. Three different 6-DOF models were proposed and their drop performances were compared with the RecurDyn model, into which MRMLG was incorporated. The dynamic models proposed in this work include the nonlinear aircraft model (NLAM), which was derived from a previous study [21]; the linear approximated aircraft model (LAAM), which linearized the nonlinear equations of motion in NLAM; and the fully approximated aircraft model (FAAM), which linearized the strut force model of the MRMLG. To validate the suitability of NLAM, LAAM, and FAAM, an MRMLG-integrated 7-DOF aircraft model was constructed using RecurDyn, which is a useful multi-body dynamics analysis tool. A comparison and error analysis with the RecurDyn model were performed using the coefficient of determination (R²) and root mean square error (RMSE). The calculation of R² values and RMSE provided a comprehensive evaluation of the model’s accuracy. Furthermore, the construction of an ensemble of time series data for various drop scenarios provided valuable insights into the performance of the derived models.

It has been shown from the evaluation that the NLAM model demonstrates a high level of accuracy, with R² values consistently above 0.93 for all variables, closely matching the results of the RecurDyn model. This indicates that NLAM is a reliable choice for simulating the system’s behavior. The LAAM model showed a satisfactory displacement response, with R² values above 0.85 and force model accuracy above 0.7, making it suitable for drop simulation. However, there were instances of data-shifting observed in the time series data, affecting the accuracy of certain variables. On the other hand, the FAAM model exhibited poor accuracy, with R² values below 0.45 for all variables except hydraulic force. This indicates that FAAM is not suitable for drop test modeling due to its unreliable performance.

The RMSE analysis further supports the following findings. NLAM consistently achieved low RMSE values, with displacements within 2 mm and angular displacements below 0.1 degrees. LAAM showed slightly higher RMSE values, but was still within acceptable ranges for practical applications. However, FAAM displayed significantly higher RMSE values, indicating large errors in its predictions. In summary, based on the ensemble analysis, NLAM has emerged as the most accurate model, closely matching the results of RecurDyn. LAAM showed a satisfactory performance, with some limitations due to data-shifting, while FAAM demonstrated poor accuracy. These findings highlight the importance of selecting an appropriate model for drop test simulations in order to ensure that reliable and accurate results are achieved for an aircraft landing gear system.